diff options
Diffstat (limited to 'boost/geometry')
88 files changed, 3857 insertions, 2243 deletions
diff --git a/boost/geometry/algorithms/assign.hpp b/boost/geometry/algorithms/assign.hpp index b17e305f2a..e3d664de38 100644 --- a/boost/geometry/algorithms/assign.hpp +++ b/boost/geometry/algorithms/assign.hpp @@ -24,7 +24,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/variant/apply_visitor.hpp> #include <boost/variant/static_visitor.hpp> diff --git a/boost/geometry/algorithms/buffer.hpp b/boost/geometry/algorithms/buffer.hpp index e776072a2f..5dfe9d8846 100644 --- a/boost/geometry/algorithms/buffer.hpp +++ b/boost/geometry/algorithms/buffer.hpp @@ -277,7 +277,7 @@ inline void buffer(GeometryIn const& geometry_in, rescale_policy_type rescale_policy = boost::geometry::get_rescale_policy<rescale_policy_type>(box); - detail::buffer::buffer_inserter<polygon_type>(geometry_in, std::back_inserter(geometry_out), + detail::buffer::buffer_inserter<polygon_type>(geometry_in, range::back_inserter(geometry_out), distance_strategy, side_strategy, join_strategy, diff --git a/boost/geometry/algorithms/convex_hull.hpp b/boost/geometry/algorithms/convex_hull.hpp index 71c18bb5f3..19d28bc7b5 100644 --- a/boost/geometry/algorithms/convex_hull.hpp +++ b/boost/geometry/algorithms/convex_hull.hpp @@ -83,7 +83,7 @@ struct hull_to_geometry geometry::point_order<OutputGeometry>::value, geometry::closure<OutputGeometry>::value >::apply(geometry, - std::back_inserter( + range::back_inserter( // Handle linestring, ring and polygon the same: detail::as_range < 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) diff --git a/boost/geometry/algorithms/difference.hpp b/boost/geometry/algorithms/difference.hpp index 780436f015..dc31c7db4d 100644 --- a/boost/geometry/algorithms/difference.hpp +++ b/boost/geometry/algorithms/difference.hpp @@ -159,7 +159,7 @@ inline void difference(Geometry1 const& geometry1, detail::difference::difference_insert<geometry_out>( geometry1, geometry2, robust_policy, - std::back_inserter(output_collection)); + range::back_inserter(output_collection)); } diff --git a/boost/geometry/algorithms/length.hpp b/boost/geometry/algorithms/length.hpp index 58324bcc80..cf5234da1a 100644 --- a/boost/geometry/algorithms/length.hpp +++ b/boost/geometry/algorithms/length.hpp @@ -24,8 +24,6 @@ #include <boost/concept_check.hpp> #include <boost/core/ignore_unused.hpp> -#include <boost/range.hpp> - #include <boost/mpl/fold.hpp> #include <boost/mpl/greater.hpp> #include <boost/mpl/if.hpp> @@ -34,8 +32,10 @@ #include <boost/mpl/set.hpp> #include <boost/mpl/size.hpp> #include <boost/mpl/transform.hpp> -#include <boost/type_traits.hpp> - +#include <boost/range/begin.hpp> +#include <boost/range/end.hpp> +#include <boost/range/iterator.hpp> +#include <boost/range/value_type.hpp> #include <boost/variant/apply_visitor.hpp> #include <boost/variant/static_visitor.hpp> #include <boost/variant/variant_fwd.hpp> diff --git a/boost/geometry/algorithms/remove_spikes.hpp b/boost/geometry/algorithms/remove_spikes.hpp index a96a2c29d0..080db92f6d 100644 --- a/boost/geometry/algorithms/remove_spikes.hpp +++ b/boost/geometry/algorithms/remove_spikes.hpp @@ -139,7 +139,7 @@ struct range_remove_spikes // Copy output geometry::clear(range); - std::copy(cleaned.begin(), cleaned.end(), std::back_inserter(range)); + std::copy(cleaned.begin(), cleaned.end(), range::back_inserter(range)); } }; diff --git a/boost/geometry/algorithms/simplify.hpp b/boost/geometry/algorithms/simplify.hpp index 9dd87f61a9..0b28eb7d15 100644 --- a/boost/geometry/algorithms/simplify.hpp +++ b/boost/geometry/algorithms/simplify.hpp @@ -77,7 +77,7 @@ struct simplify_copy { std::copy ( - boost::begin(range), boost::end(range), std::back_inserter(out) + boost::begin(range), boost::end(range), geometry::range::back_inserter(out) ); } }; @@ -113,7 +113,7 @@ struct simplify_range { simplify_range_insert::apply ( - range, std::back_inserter(out), max_distance, strategy + range, geometry::range::back_inserter(out), max_distance, strategy ); } } diff --git a/boost/geometry/algorithms/sym_difference.hpp b/boost/geometry/algorithms/sym_difference.hpp index 6c49ca634c..77a20c0885 100644 --- a/boost/geometry/algorithms/sym_difference.hpp +++ b/boost/geometry/algorithms/sym_difference.hpp @@ -332,7 +332,7 @@ inline void sym_difference(Geometry1 const& geometry1, detail::sym_difference::sym_difference_insert<geometry_out>( geometry1, geometry2, robust_policy, - std::back_inserter(output_collection)); + range::back_inserter(output_collection)); } diff --git a/boost/geometry/algorithms/transform.hpp b/boost/geometry/algorithms/transform.hpp index c7f3bd6cbb..f6748b11e3 100644 --- a/boost/geometry/algorithms/transform.hpp +++ b/boost/geometry/algorithms/transform.hpp @@ -162,7 +162,7 @@ struct transform_polygon geometry::clear(poly2); if (!transform_range_out<point2_type>(geometry::exterior_ring(poly1), - std::back_inserter(geometry::exterior_ring(poly2)), strategy)) + range::back_inserter(geometry::exterior_ring(poly2)), strategy)) { return false; } @@ -189,7 +189,7 @@ struct transform_polygon for ( ; it1 != boost::end(rings1); ++it1, ++it2) { if ( ! transform_range_out<point2_type>(*it1, - std::back_inserter(*it2), + range::back_inserter(*it2), strategy) ) { return false; @@ -228,7 +228,7 @@ struct transform_range // Should NOT be done here! // geometry::clear(range2); return transform_range_out<point_type>(range1, - std::back_inserter(range2), strategy); + range::back_inserter(range2), strategy); } }; diff --git a/boost/geometry/algorithms/union.hpp b/boost/geometry/algorithms/union.hpp index ff16c60ea1..d15ee7655f 100644 --- a/boost/geometry/algorithms/union.hpp +++ b/boost/geometry/algorithms/union.hpp @@ -271,7 +271,7 @@ inline void union_(Geometry1 const& geometry1, concept::check<geometry_out>(); detail::union_::union_insert<geometry_out>(geometry1, geometry2, - std::back_inserter(output_collection)); + range::back_inserter(output_collection)); } diff --git a/boost/geometry/core/access.hpp b/boost/geometry/core/access.hpp index 57f925763d..866fb19a4d 100644 --- a/boost/geometry/core/access.hpp +++ b/boost/geometry/core/access.hpp @@ -17,9 +17,10 @@ #include <cstddef> +#include <boost/core/ignore_unused.hpp> #include <boost/mpl/assert.hpp> -#include <boost/concept_check.hpp> #include <boost/type_traits/is_pointer.hpp> +#include <boost/type_traits/remove_pointer.hpp> #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/point_type.hpp> @@ -271,7 +272,7 @@ inline typename coordinate_type<Geometry>::type get(Geometry const& geometry #endif ) { - boost::ignore_unused_variable_warning(dummy); + boost::ignore_unused(dummy); typedef core_dispatch::access < @@ -306,7 +307,7 @@ inline void set(Geometry& geometry #endif ) { - boost::ignore_unused_variable_warning(dummy); + boost::ignore_unused(dummy); typedef core_dispatch::access < @@ -341,7 +342,7 @@ inline typename coordinate_type<Geometry>::type get(Geometry const& geometry #endif ) { - boost::ignore_unused_variable_warning(dummy); + boost::ignore_unused(dummy); typedef core_dispatch::indexed_access < @@ -378,7 +379,7 @@ inline void set(Geometry& geometry #endif ) { - boost::ignore_unused_variable_warning(dummy); + boost::ignore_unused(dummy); typedef core_dispatch::indexed_access < diff --git a/boost/geometry/core/closure.hpp b/boost/geometry/core/closure.hpp index b69dcda140..405a5a5f35 100644 --- a/boost/geometry/core/closure.hpp +++ b/boost/geometry/core/closure.hpp @@ -22,11 +22,11 @@ #include <boost/mpl/assert.hpp> #include <boost/mpl/size_t.hpp> #include <boost/range/value_type.hpp> -#include <boost/type_traits/remove_const.hpp> #include <boost/geometry/core/ring_type.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry { diff --git a/boost/geometry/core/coordinate_dimension.hpp b/boost/geometry/core/coordinate_dimension.hpp index 85da6b424e..19022f2fd2 100644 --- a/boost/geometry/core/coordinate_dimension.hpp +++ b/boost/geometry/core/coordinate_dimension.hpp @@ -18,7 +18,6 @@ #include <cstddef> #include <boost/mpl/assert.hpp> -#include <boost/mpl/equal_to.hpp> #include <boost/static_assert.hpp> #include <boost/geometry/core/point_type.hpp> @@ -94,13 +93,7 @@ struct dimension template <typename Geometry, int Dimensions> inline void assert_dimension() { - BOOST_STATIC_ASSERT(( - boost::mpl::equal_to - < - boost::mpl::int_<geometry::dimension<Geometry>::value>, - boost::mpl::int_<Dimensions> - >::type::value - )); + BOOST_STATIC_ASSERT(( static_cast<int>(dimension<Geometry>::value) == Dimensions )); } /*! diff --git a/boost/geometry/core/coordinate_type.hpp b/boost/geometry/core/coordinate_type.hpp index f138d59c8a..47c028fec2 100644 --- a/boost/geometry/core/coordinate_type.hpp +++ b/boost/geometry/core/coordinate_type.hpp @@ -17,8 +17,8 @@ #include <boost/mpl/assert.hpp> -#include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/tag.hpp> #include <boost/geometry/util/bare_type.hpp> #include <boost/geometry/util/promote_floating_point.hpp> diff --git a/boost/geometry/core/exterior_ring.hpp b/boost/geometry/core/exterior_ring.hpp index e5c97dd301..fa5a0a2b98 100644 --- a/boost/geometry/core/exterior_ring.hpp +++ b/boost/geometry/core/exterior_ring.hpp @@ -17,6 +17,7 @@ #include <boost/mpl/assert.hpp> +#include <boost/type_traits/is_const.hpp> #include <boost/type_traits/remove_const.hpp> diff --git a/boost/geometry/core/interior_type.hpp b/boost/geometry/core/interior_type.hpp index e6c4537e41..66544b79af 100644 --- a/boost/geometry/core/interior_type.hpp +++ b/boost/geometry/core/interior_type.hpp @@ -17,7 +17,10 @@ #include <boost/mpl/assert.hpp> +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_const.hpp> #include <boost/type_traits/remove_const.hpp> +#include <boost/type_traits/remove_reference.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> diff --git a/boost/geometry/core/mutable_range.hpp b/boost/geometry/core/mutable_range.hpp index 9b53e40577..7cb5137830 100644 --- a/boost/geometry/core/mutable_range.hpp +++ b/boost/geometry/core/mutable_range.hpp @@ -17,8 +17,8 @@ #include <cstddef> +#include <boost/range/value_type.hpp> #include <boost/type_traits/remove_reference.hpp> -#include <boost/range.hpp> namespace boost { namespace geometry diff --git a/boost/geometry/core/point_order.hpp b/boost/geometry/core/point_order.hpp index d43adbd03a..df6d4bbfff 100644 --- a/boost/geometry/core/point_order.hpp +++ b/boost/geometry/core/point_order.hpp @@ -22,11 +22,11 @@ #include <boost/mpl/assert.hpp> #include <boost/range/value_type.hpp> -#include <boost/type_traits/remove_const.hpp> #include <boost/geometry/core/ring_type.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry { diff --git a/boost/geometry/core/radius.hpp b/boost/geometry/core/radius.hpp index ee6d5d246f..6cbf6ac8f2 100644 --- a/boost/geometry/core/radius.hpp +++ b/boost/geometry/core/radius.hpp @@ -24,6 +24,7 @@ #include <cstddef> #include <boost/static_assert.hpp> +#include <boost/type_traits/is_pointer.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> diff --git a/boost/geometry/core/reverse_dispatch.hpp b/boost/geometry/core/reverse_dispatch.hpp index 2e4fb8005f..c95182a34d 100644 --- a/boost/geometry/core/reverse_dispatch.hpp +++ b/boost/geometry/core/reverse_dispatch.hpp @@ -18,9 +18,8 @@ #include <cstddef> -#include <boost/type_traits.hpp> #include <boost/mpl/if.hpp> -#include <boost/mpl/greater.hpp> +#include <boost/type_traits/integral_constant.hpp> #include <boost/geometry/core/geometry_id.hpp> diff --git a/boost/geometry/core/tag.hpp b/boost/geometry/core/tag.hpp index 6bffcb932d..14a2242708 100644 --- a/boost/geometry/core/tag.hpp +++ b/boost/geometry/core/tag.hpp @@ -14,7 +14,6 @@ #ifndef BOOST_GEOMETRY_CORE_TAG_HPP #define BOOST_GEOMETRY_CORE_TAG_HPP -#include <boost/mpl/assert.hpp> #include <boost/geometry/core/tags.hpp> #include <boost/geometry/util/bare_type.hpp> diff --git a/boost/geometry/core/tag_cast.hpp b/boost/geometry/core/tag_cast.hpp index cafd13cba9..e03d841f49 100644 --- a/boost/geometry/core/tag_cast.hpp +++ b/boost/geometry/core/tag_cast.hpp @@ -16,7 +16,7 @@ #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_base_of.hpp> namespace boost { namespace geometry { diff --git a/boost/geometry/extensions/algorithms/inverse.hpp b/boost/geometry/extensions/algorithms/inverse.hpp new file mode 100644 index 0000000000..44f415be42 --- /dev/null +++ b/boost/geometry/extensions/algorithms/inverse.hpp @@ -0,0 +1,59 @@ +// 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_INVERSE_HPP +#define BOOST_GEOMETRY_ALGORITHMS_INVERSE_HPP + +#include <boost/geometry.hpp> + +namespace boost { namespace geometry +{ + +// TODO: this is shared with sectionalize, move to somewhere else (assign?) +template <typename Box, typename Value> +inline void enlarge_box(Box& box, Value value) +{ + geometry::set<0, 0>(box, geometry::get<0, 0>(box) - value); + geometry::set<0, 1>(box, geometry::get<0, 1>(box) - value); + geometry::set<1, 0>(box, geometry::get<1, 0>(box) + value); + geometry::set<1, 1>(box, geometry::get<1, 1>(box) + value); +} + +// TODO: when this might be moved outside extensions it should of course +// input/output a Geometry, instead of a WKT +template <typename Geometry, typename Value> +inline std::string inverse(std::string const& wkt, Value margin) +{ + Geometry geometry; + read_wkt(wkt, geometry); + + geometry::correct(geometry); + + geometry::model::box<typename point_type<Geometry>::type> env; + geometry::envelope(geometry, env); + + // Make its envelope a bit larger + enlarge_box(env, margin); + + Geometry env_as_polygon; + geometry::convert(env, env_as_polygon); + + Geometry inversed_result; + geometry::difference(env_as_polygon, geometry, inversed_result); + + std::ostringstream out; + + out << geometry::wkt(inversed_result); + return out.str(); +} + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_ALGORITHMS_INVERSE_HPP diff --git a/boost/geometry/geometries/adapted/boost_polygon/ring_proxy.hpp b/boost/geometry/geometries/adapted/boost_polygon/ring_proxy.hpp index 1616369da1..84f72aadcb 100644 --- a/boost/geometry/geometries/adapted/boost_polygon/ring_proxy.hpp +++ b/boost/geometry/geometries/adapted/boost_polygon/ring_proxy.hpp @@ -55,7 +55,7 @@ template <> struct modify<false> { template <typename Ring, typename Point> - static inline void push_back(Ring& ring, Point const& point) + static inline void push_back(Ring& /*ring*/, Point const& /*point*/) { } diff --git a/boost/geometry/index/detail/rtree/pack_create.hpp b/boost/geometry/index/detail/rtree/pack_create.hpp index e56ce076d2..d1491b8d47 100644 --- a/boost/geometry/index/detail/rtree/pack_create.hpp +++ b/boost/geometry/index/detail/rtree/pack_create.hpp @@ -200,6 +200,13 @@ private: {} template <typename Indexable> + explicit expandable_box(Indexable const& indexable) + : m_initialized(true) + { + detail::bounds(indexable, m_box); + } + + template <typename Indexable> void expand(Indexable const& indexable) { if ( !m_initialized ) @@ -260,9 +267,12 @@ private: // reserve space for values rtree::elements(l).reserve(values_count); // MAY THROW (A) + // calculate values box and copy values - expandable_box<Box> elements_box; - for ( ; first != last ; ++first ) + // initialize the box explicitly to avoid GCC-4.4 uninitialized variable warnings with O2 + expandable_box<Box> elements_box(translator(*(first->second))); + rtree::elements(l).push_back(*(first->second)); // MAY THROW (A?,C) + for ( ++first ; first != last ; ++first ) { // NOTE: push_back() must be called at the end in order to support move_iterator. // The iterator is dereferenced 2x (no temporary reference) to support diff --git a/boost/geometry/index/rtree.hpp b/boost/geometry/index/rtree.hpp index ea7fc74ed3..6d3704ff8c 100644 --- a/boost/geometry/index/rtree.hpp +++ b/boost/geometry/index/rtree.hpp @@ -1263,16 +1263,16 @@ public: inline bounds_type bounds() const { bounds_type result; - if ( !m_members.root ) + // in order to suppress the uninitialized variable warnings + geometry::assign_inverse(result); + + if ( m_members.root ) { - geometry::assign_inverse(result); - return result; + detail::rtree::visitors::children_box<value_type, options_type, translator_type, box_type, allocators_type> + box_v(result, m_members.translator()); + detail::rtree::apply_visitor(box_v, *m_members.root); } - detail::rtree::visitors::children_box<value_type, options_type, translator_type, box_type, allocators_type> - box_v(result, m_members.translator()); - detail::rtree::apply_visitor(box_v, *m_members.root); - return result; } diff --git a/boost/geometry/io/svg/svg_mapper.hpp b/boost/geometry/io/svg/svg_mapper.hpp index c8e63d5ab7..e06f2acc24 100644 --- a/boost/geometry/io/svg/svg_mapper.hpp +++ b/boost/geometry/io/svg/svg_mapper.hpp @@ -2,10 +2,11 @@ // Copyright (c) 2009-2015 Barend Gehrels, Amsterdam, the Netherlands. -// 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 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -90,26 +91,36 @@ struct svg_map<point_tag, Point> } }; -template <typename Box> -struct svg_map<box_tag, Box> +template <typename BoxSeg1, typename BoxSeg2> +struct svg_map_box_seg { template <typename TransformStrategy> static inline void apply(std::ostream& stream, std::string const& style, int size, - Box const& box, TransformStrategy const& strategy) + BoxSeg1 const& box_seg, TransformStrategy const& strategy) { - model::box<detail::svg::svg_point_type> ibox; + BoxSeg2 ibox_seg; - // Fix bug in gcc compiler warning for possible uninitialation + // Fix bug in gcc compiler warning for possible uninitialization #if defined(BOOST_GCC) - geometry::assign_zero(ibox); + geometry::assign_zero(ibox_seg); #endif - geometry::transform(box, ibox, strategy); + geometry::transform(box_seg, ibox_seg, strategy); - stream << geometry::svg(ibox, style, size) << std::endl; + stream << geometry::svg(ibox_seg, style, size) << std::endl; } }; +template <typename Box> +struct svg_map<box_tag, Box> + : svg_map_box_seg<Box, model::box<detail::svg::svg_point_type> > +{}; + +template <typename Segment> +struct svg_map<segment_tag, Segment> + : svg_map_box_seg<Segment, model::segment<detail::svg::svg_point_type> > +{}; + template <typename Range1, typename Range2> struct svg_map_range @@ -125,25 +136,6 @@ struct svg_map_range } }; -template <typename Segment> -struct svg_map<segment_tag, Segment> -{ - template <typename TransformStrategy> - static inline void apply(std::ostream& stream, - std::string const& style, int size, - Segment const& segment, TransformStrategy const& strategy) - { - typedef segment_view<Segment> view_type; - view_type range(segment); - svg_map_range - < - view_type, - model::linestring<detail::svg::svg_point_type> - >::apply(stream, style, size, range, strategy); - } -}; - - template <typename Ring> struct svg_map<ring_tag, Ring> : svg_map_range<Ring, model::ring<detail::svg::svg_point_type> > diff --git a/boost/geometry/io/svg/write_svg.hpp b/boost/geometry/io/svg/write_svg.hpp index fba3cbebaf..4a518a0815 100644 --- a/boost/geometry/io/svg/write_svg.hpp +++ b/boost/geometry/io/svg/write_svg.hpp @@ -3,6 +3,11 @@ // Copyright (c) 2009-2012 Barend Gehrels, Amsterdam, the Netherlands. // Copyright (c) 2014 Adam Wulkiewicz, Lodz, Poland. +// 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 + // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -66,9 +71,9 @@ struct svg_box typedef typename coordinate_type<Box>::type ct; ct x = geometry::get<geometry::min_corner, 0>(box); ct y = geometry::get<geometry::min_corner, 1>(box); - ct width = max BOOST_PREVENT_MACRO_SUBSTITUTION(1, + ct width = max BOOST_PREVENT_MACRO_SUBSTITUTION (ct(1), geometry::get<geometry::max_corner, 0>(box) - x); - ct height = max BOOST_PREVENT_MACRO_SUBSTITUTION (1, + ct height = max BOOST_PREVENT_MACRO_SUBSTITUTION (ct(1), geometry::get<geometry::max_corner, 1>(box) - y); os << "<rect x=\"" << x << "\" y=\"" << y @@ -77,6 +82,24 @@ struct svg_box } }; +template <typename Segment> +struct svg_segment +{ + template <typename Char, typename Traits> + static inline void apply(std::basic_ostream<Char, Traits>& os, + Segment const& segment, std::string const& style, int) + { + typedef typename coordinate_type<Segment>::type ct; + ct x1 = geometry::get<0, 0>(segment); + ct y1 = geometry::get<0, 1>(segment); + ct x2 = geometry::get<1, 0>(segment); + ct y2 = geometry::get<1, 1>(segment); + + os << "<line x1=\"" << x1 << "\" y1=\"" << y1 + << "\" x2=\"" << x2 << "\" y2=\"" << y2 + << "\" style=\"" << style << "\"/>"; + } +}; /*! \brief Stream ranges as SVG @@ -204,6 +227,9 @@ struct svg template <typename Point> struct svg<point_tag, Point> : detail::svg::svg_point<Point> {}; +template <typename Segment> +struct svg<segment_tag, Segment> : detail::svg::svg_segment<Segment> {}; + template <typename Box> struct svg<box_tag, Box> : detail::svg::svg_box<Box> {}; diff --git a/boost/geometry/io/svg/write_svg_multi.hpp b/boost/geometry/io/svg/write_svg_multi.hpp index a794120c06..ff349ebd75 100644 --- a/boost/geometry/io/svg/write_svg_multi.hpp +++ b/boost/geometry/io/svg/write_svg_multi.hpp @@ -2,6 +2,11 @@ // Copyright (c) 2009-2012 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 + // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -53,6 +58,32 @@ struct svg_multi namespace dispatch { +template <typename MultiPoint> +struct svg<multi_point_tag, MultiPoint> + : detail::svg::svg_multi + < + MultiPoint, + detail::svg::svg_point + < + typename boost::range_value<MultiPoint>::type + > + + > +{}; + +template <typename MultiLinestring> +struct svg<multi_linestring_tag, MultiLinestring> + : detail::svg::svg_multi + < + MultiLinestring, + detail::svg::svg_range + < + typename boost::range_value<MultiLinestring>::type, + detail::svg::prefix_linestring + > + + > +{}; template <typename MultiPolygon> struct svg<multi_polygon_tag, MultiPolygon> diff --git a/boost/geometry/io/wkt/read.hpp b/boost/geometry/io/wkt/read.hpp index a39c8f89d4..7924b70283 100644 --- a/boost/geometry/io/wkt/read.hpp +++ b/boost/geometry/io/wkt/read.hpp @@ -27,9 +27,12 @@ #include <boost/algorithm/string.hpp> #include <boost/mpl/if.hpp> -#include <boost/range.hpp> - -#include <boost/type_traits.hpp> +#include <boost/range/begin.hpp> +#include <boost/range/end.hpp> +#include <boost/range/size.hpp> +#include <boost/range/value_type.hpp> +#include <boost/type_traits/is_same.hpp> +#include <boost/type_traits/remove_reference.hpp> #include <boost/geometry/algorithms/assign.hpp> #include <boost/geometry/algorithms/append.hpp> diff --git a/boost/geometry/policies/robustness/get_rescale_policy.hpp b/boost/geometry/policies/robustness/get_rescale_policy.hpp index e7bef1d2e1..daf32d3057 100644 --- a/boost/geometry/policies/robustness/get_rescale_policy.hpp +++ b/boost/geometry/policies/robustness/get_rescale_policy.hpp @@ -20,8 +20,9 @@ #include <cstddef> -#include <boost/type_traits.hpp> #include <boost/mpl/assert.hpp> +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/type_traits/is_same.hpp> #include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/tag_cast.hpp> @@ -252,9 +253,15 @@ struct rescale_policy_type false #else boost::is_floating_point - < - typename geometry::coordinate_type<Point>::type - >::type::value + < + typename geometry::coordinate_type<Point>::type + >::type::value + && + boost::is_same + < + typename geometry::coordinate_system<Point>::type, + geometry::cs::cartesian + >::value #endif > { diff --git a/boost/geometry/policies/robustness/rescale_policy.hpp b/boost/geometry/policies/robustness/rescale_policy.hpp index 8f34bf3dee..b92f6e1ec7 100644 --- a/boost/geometry/policies/robustness/rescale_policy.hpp +++ b/boost/geometry/policies/robustness/rescale_policy.hpp @@ -18,8 +18,6 @@ #include <cstddef> -#include <boost/type_traits.hpp> - #include <boost/geometry/policies/robustness/segment_ratio.hpp> #include <boost/geometry/policies/robustness/segment_ratio_type.hpp> #include <boost/geometry/policies/robustness/robust_point_type.hpp> diff --git a/boost/geometry/policies/robustness/robust_type.hpp b/boost/geometry/policies/robustness/robust_type.hpp index 4cfb2c4d91..7342f905bf 100644 --- a/boost/geometry/policies/robustness/robust_type.hpp +++ b/boost/geometry/policies/robustness/robust_type.hpp @@ -13,8 +13,8 @@ #define BOOST_GEOMETRY_POLICIES_ROBUSTNESS_ROBUST_TYPE_HPP -#include <boost/type_traits.hpp> #include <boost/config.hpp> +#include <boost/type_traits/is_floating_point.hpp> namespace boost { namespace geometry diff --git a/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp b/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp index 968f9fecb9..6740b519e2 100644 --- a/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp +++ b/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp @@ -55,6 +55,10 @@ struct decide_covered_by }; +// WARNING +// This strategy is not suitable for boxes in non-cartesian CSes having edges +// longer than 180deg because e.g. the SSF formula picks the side of the closer +// longitude, so for long edges the side is the opposite. template <typename Point, typename Box, typename Decide = decide_within> struct point_in_box_by_side { @@ -93,58 +97,6 @@ struct point_in_box_by_side } // namespace within -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - - -namespace within { namespace services -{ - -template <typename Point, typename Box> -struct default_strategy - < - point_tag, box_tag, - point_tag, areal_tag, - spherical_tag, spherical_tag, - Point, Box - > -{ - typedef within::point_in_box_by_side - < - Point, Box, within::decide_within - > type; -}; - - - -}} // namespace within::services - - -namespace covered_by { namespace services -{ - - -template <typename Point, typename Box> -struct default_strategy - < - point_tag, box_tag, - point_tag, areal_tag, - spherical_tag, spherical_tag, - Point, Box - > -{ - typedef within::point_in_box_by_side - < - Point, Box, within::decide_covered_by - > type; -}; - - -}} // namespace covered_by::services - - -#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - - }}} // namespace boost::geometry::strategy diff --git a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp index 641533fc6a..9e2ec2c4ff 100644 --- a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp @@ -3,8 +3,9 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2013, 2014. -// Modifications copyright (c) 2013, 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2013, 2014, 2016. +// Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -13,8 +14,6 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - #ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP @@ -35,15 +34,34 @@ namespace boost { namespace geometry namespace strategy { namespace within { +// 1 deg or pi/180 rad +template <typename Point, + typename CalculationType = typename coordinate_type<Point>::type> +struct winding_small_angle +{ + typedef typename coordinate_system<Point>::type cs_t; + typedef math::detail::constants_on_spheroid + < + CalculationType, + typename cs_t::units + > constants; + + static inline CalculationType apply() + { + return constants::half_period() / CalculationType(180); + } +}; + // Fix for https://svn.boost.org/trac/boost/ticket/9628 -// For floating point coordinates, the <1> coordinate of a point is compared +// For floating point coordinates, the <D> coordinate of a point is compared // with the segment's points using some EPS. If the coordinates are "equal" // the sides are calculated. Therefore we can treat a segment as a long areal // geometry having some width. There is a small ~triangular area somewhere // between the segment's effective area and a segment's line used in sides // calculation where the segment is on the one side of the line but on the // other side of a segment (due to the width). +// Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP. // For the s1 of a segment going NE the real side is RIGHT but the point may // be detected as LEFT, like this: // RIGHT @@ -54,6 +72,9 @@ namespace strategy { namespace within // _____7 // _____/ // _____/ +// In the code below actually D = 0, so segments are nearly-vertical +// Called when the point is on the same level as one of the segment's points +// but the point is not aligned with a vertical segment template <typename CSTag> struct winding_side_equal { @@ -62,122 +83,228 @@ struct winding_side_equal CSTag >::type strategy_side_type; - template <size_t D, typename Point, typename PointOfSegment> + template <typename Point, typename PointOfSegment> static inline int apply(Point const& point, PointOfSegment const& se, int count) { - // Create a vertical segment intersecting the original segment's endpoint - // equal to the point, with the derived direction (UP/DOWN). - // Set only the 2 first coordinates, the other ones are ignored + typedef typename coordinate_type<PointOfSegment>::type scoord_t; + typedef typename coordinate_system<PointOfSegment>::type::units units_t; + + if (math::equals(get<1>(point), get<1>(se))) + return 0; + + // Create a horizontal segment intersecting the original segment's endpoint + // equal to the point, with the derived direction (E/W). PointOfSegment ss1, ss2; - set<1-D>(ss1, get<1-D>(se)); - set<1-D>(ss2, get<1-D>(se)); - if (count > 0) // UP + set<1>(ss1, get<1>(se)); + set<0>(ss1, get<0>(se)); + set<1>(ss2, get<1>(se)); + scoord_t ss20 = get<0>(se); + if (count > 0) { - set<D>(ss1, 0); - set<D>(ss2, 1); + ss20 += winding_small_angle<PointOfSegment>::apply(); } - else // DOWN + else { - set<D>(ss1, 1); - set<D>(ss2, 0); + ss20 -= winding_small_angle<PointOfSegment>::apply(); } + math::normalize_longitude<units_t>(ss20); + set<0>(ss2, ss20); + // Check the side using this vertical segment return strategy_side_type::apply(ss1, ss2, point); } }; - // The optimization for cartesian template <> struct winding_side_equal<cartesian_tag> { - template <size_t D, typename Point, typename PointOfSegment> + template <typename Point, typename PointOfSegment> static inline int apply(Point const& point, PointOfSegment const& se, int count) { - return math::equals(get<1-D>(point), get<1-D>(se)) ? + // NOTE: for D=0 the signs would be reversed + return math::equals(get<1>(point), get<1>(se)) ? 0 : - get<1-D>(point) < get<1-D>(se) ? + get<1>(point) < get<1>(se) ? // assuming count is equal to 1 or -1 - count : // ( count > 0 ? 1 : -1) : - -count; // ( count > 0 ? -1 : 1) ; + -count : // ( count > 0 ? -1 : 1) : + count; // ( count > 0 ? 1 : -1) ; } }; -template <typename CSTag> -struct winding_side_between +template <typename Point, + typename CalculationType, + typename CSTag = typename cs_tag<Point>::type> +struct winding_check_touch { - typedef typename strategy::side::services::default_strategy - < - CSTag - >::type strategy_side_type; + typedef CalculationType calc_t; + typedef typename coordinate_system<Point>::type::units units_t; + typedef math::detail::constants_on_spheroid<CalculationType, units_t> constants; - template <size_t D, typename Point, typename PointOfSegment> + template <typename PointOfSegment, typename State> static inline int apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - int count) + PointOfSegment const& seg1, + PointOfSegment const& seg2, + State& state, + bool& eq1, + bool& eq2) { - // Create a vertical segment intersecting the original segment's endpoint - // equal to the point, with the derived direction (UP/DOWN). - // Set only the 2 first coordinates, the other ones are ignored - PointOfSegment ss1, ss2; - set<1-D>(ss1, get<1-D>(s1)); - set<1-D>(ss2, get<1-D>(s1)); - - if (count > 0) // UP + calc_t const pi = constants::half_period(); + calc_t const pi2 = pi / calc_t(2); + + calc_t const px = get<0>(point); + calc_t const s1x = get<0>(seg1); + calc_t const s2x = get<0>(seg2); + calc_t const py = get<1>(point); + calc_t const s1y = get<1>(seg1); + calc_t const s2y = get<1>(seg2); + + // NOTE: lat in {-90, 90} and arbitrary lon + // it doesn't matter what lon it is if it's a pole + // so e.g. if one of the segment endpoints is a pole + // then only the other lon matters + + bool eq1_strict = math::equals(s1x, px); + bool eq2_strict = math::equals(s2x, px); + + eq1 = eq1_strict // lon strictly equal to s1 + || math::equals(s1y, pi2) || math::equals(s1y, -pi2); // s1 is pole + eq2 = eq2_strict // lon strictly equal to s2 + || math::equals(s2y, pi2) || math::equals(s2y, -pi2); // s2 is pole + + // segment overlapping pole + calc_t s1x_anti = s1x + constants::half_period(); + math::normalize_longitude<units_t, calc_t>(s1x_anti); + bool antipodal = math::equals(s2x, s1x_anti); + if (antipodal) { - set<D>(ss1, 0); - set<D>(ss2, 1); + eq1 = eq2 = eq1 || eq2; + + // segment overlapping pole and point is pole + if (math::equals(py, pi2) || math::equals(py, -pi2)) + { + eq1 = eq2 = true; + } } - else // DOWN + + // Both equal p -> segment vertical + // The only thing which has to be done is check if point is ON segment + if (eq1 && eq2) { - set<D>(ss1, 1); - set<D>(ss2, 0); + // segment endpoints on the same sides of the globe + if (! antipodal + // p's lat between segment endpoints' lats + ? (s1y <= py && s2y >= py) || (s2y <= py && s1y >= py) + // going through north or south pole? + : (pi - s1y - s2y <= pi + ? (eq1_strict && s1y <= py) || (eq2_strict && s2y <= py) // north + || math::equals(py, pi2) // point on north pole + : (eq1_strict && s1y >= py) || (eq2_strict && s2y >= py)) // south + || math::equals(py, -pi2) // point on south pole + ) + { + state.m_touches = true; + } + return true; } + return false; + } +}; +// The optimization for cartesian +template <typename Point, typename CalculationType> +struct winding_check_touch<Point, CalculationType, cartesian_tag> +{ + typedef CalculationType calc_t; - int const seg_side = strategy_side_type::apply(ss1, ss2, s2); + template <typename PointOfSegment, typename State> + static inline bool apply(Point const& point, + PointOfSegment const& seg1, + PointOfSegment const& seg2, + State& state, + bool& eq1, + bool& eq2) + { + calc_t const px = get<0>(point); + calc_t const s1x = get<0>(seg1); + calc_t const s2x = get<0>(seg2); - if (seg_side != 0) // segment not vertical + eq1 = math::equals(s1x, px); + eq2 = math::equals(s2x, px); + + // Both equal p -> segment vertical + // The only thing which has to be done is check if point is ON segment + if (eq1 && eq2) { - if (strategy_side_type::apply(ss1, ss2, point) == -seg_side) // point on the opposite side than s2 - { - return -seg_side; - } - else + calc_t const py = get<1>(point); + calc_t const s1y = get<1>(seg1); + calc_t const s2y = get<1>(seg2); + if ((s1y <= py && s2y >= py) || (s2y <= py && s1y >= py)) { - set<1-D>(ss1, get<1-D>(s2)); - set<1-D>(ss2, get<1-D>(s2)); - - if (strategy_side_type::apply(ss1, ss2, point) == seg_side) // point behind s2 - { - return seg_side; - } + state.m_touches = true; } + return true; } - - // segment is vertical or point is between p1 and p2 - return strategy_side_type::apply(s1, s2, point); + return false; } }; -// The specialization for cartesian -template <> -struct winding_side_between<cartesian_tag> + +// Called if point is not aligned with a vertical segment +template <typename Point, + typename CalculationType, + typename CSTag = typename cs_tag<Point>::type> +struct winding_calculate_count { - typedef strategy::side::services::default_strategy - < - cartesian_tag - >::type strategy_side_type; + typedef CalculationType calc_t; + typedef typename coordinate_system<Point>::type::units units_t; - template <size_t D, typename Point, typename PointOfSegment> - static inline int apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - int /*count*/) + static inline bool greater(calc_t const& l, calc_t const& r) + { + calc_t diff = l - r; + math::normalize_longitude<units_t, calc_t>(diff); + return diff > calc_t(0); + } + + static inline int apply(calc_t const& p, + calc_t const& s1, calc_t const& s2, + bool eq1, bool eq2) + { + // Probably could be optimized by avoiding normalization for some comparisons + // e.g. s1 > p could be calculated from p > s1 + + // If both segment endpoints were poles below checks wouldn't be enough + // but this means that either both are the same or that they are N/S poles + // and therefore the segment is not valid. + // If needed (eq1 && eq2 ? 0) could be returned + + return + eq1 ? (greater(s2, p) ? 1 : -1) // Point on level s1, E/W depending on s2 + : eq2 ? (greater(s1, p) ? -1 : 1) // idem + : greater(p, s1) && greater(s2, p) ? 2 // Point between s1 -> s2 --> E + : greater(p, s2) && greater(s1, p) ? -2 // Point between s2 -> s1 --> W + : 0; + } +}; +// The optimization for cartesian +template <typename Point, typename CalculationType> +struct winding_calculate_count<Point, CalculationType, cartesian_tag> +{ + typedef CalculationType calc_t; + + static inline int apply(calc_t const& p, + calc_t const& s1, calc_t const& s2, + bool eq1, bool eq2) { - return strategy_side_type::apply(s1, s2, point); + return + eq1 ? (s2 > p ? 1 : -1) // Point on level s1, E/W depending on s2 + : eq2 ? (s1 > p ? -1 : 1) // idem + : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> E + : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> W + : 0; } }; @@ -234,6 +361,9 @@ class winding public : friend class winding; + template <typename P, typename CT, typename CST> + friend struct winding_check_touch; + inline counter() : m_count(0) , m_touches(false) @@ -242,48 +372,21 @@ class winding }; - template <size_t D> - static inline int check_touch(Point const& point, - PointOfSegment const& seg1, PointOfSegment const& seg2, - counter& state) - { - calculation_type const p = get<D>(point); - calculation_type const s1 = get<D>(seg1); - calculation_type const s2 = get<D>(seg2); - if ((s1 <= p && s2 >= p) || (s2 <= p && s1 >= p)) - { - state.m_touches = true; - } - return 0; - } - - - template <size_t D> static inline int check_segment(Point const& point, PointOfSegment const& seg1, PointOfSegment const& seg2, counter& state, bool& eq1, bool& eq2) { - calculation_type const p = get<D>(point); - calculation_type const s1 = get<D>(seg1); - calculation_type const s2 = get<D>(seg2); - - // Check if one of segment endpoints is at same level of point - eq1 = math::equals(s1, p); - eq2 = math::equals(s2, p); - - if (eq1 && eq2) + if (winding_check_touch<Point, calculation_type> + ::apply(point, seg1, seg2, state, eq1, eq2)) { - // Both equal p -> segment is horizontal (or vertical for D=0) - // The only thing which has to be done is check if point is ON segment - return check_touch<1 - D>(point, seg1, seg2, state); + return 0; } - return - eq1 ? (s2 > p ? 1 : -1) // Point on level s1, UP/DOWN depending on s2 - : eq2 ? (s1 > p ? -1 : 1) // idem - : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> UP - : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> DOWN - : 0; + calculation_type const p = get<0>(point); + calculation_type const s1 = get<0>(seg1); + calculation_type const s2 = get<0>(seg2); + return winding_calculate_count<Point, calculation_type> + ::apply(p, s1, s2, eq1, eq2); } @@ -304,19 +407,18 @@ public : bool eq2 = false; boost::ignore_unused(eq2); - int count = check_segment<1>(point, s1, s2, state, eq1, eq2); + int count = check_segment(point, s1, s2, state, eq1, eq2); if (count != 0) { int side = 0; if (count == 1 || count == -1) { - side = winding_side_equal<cs_t> - ::template apply<1>(point, eq1 ? s1 : s2, count); + side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count); } else // count == 2 || count == -2 { - side = winding_side_between<cs_t> - ::template apply<1>(point, s1, s2, count); + // 1 left, -1 right + side = strategy_side_type::apply(s1, s2, point); } if (side == 0) diff --git a/boost/geometry/strategies/cartesian/area_surveyor.hpp b/boost/geometry/strategies/cartesian/area_surveyor.hpp index 74b63532c0..ba76f67946 100644 --- a/boost/geometry/strategies/cartesian/area_surveyor.hpp +++ b/boost/geometry/strategies/cartesian/area_surveyor.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. +// 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 + // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -17,7 +22,7 @@ #include <boost/mpl/if.hpp> -#include <boost/geometry/arithmetic/determinant.hpp> +//#include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/coordinate_dimension.hpp> #include <boost/geometry/util/select_most_precise.hpp> @@ -98,8 +103,20 @@ public : PointOfSegment const& p2, summation& state) { + // Below formulas are equivalent, however the two lower ones + // suffer less from accuracy loss for great values of coordinates. + // See: https://svn.boost.org/trac/boost/ticket/11928 + // SUM += x2 * y1 - x1 * y2; - state.sum += detail::determinant<return_type>(p2, p1); + // state.sum += detail::determinant<return_type>(p2, p1); + + // SUM += (x2 - x1) * (y2 + y1) + //state.sum += (return_type(get<0>(p2)) - return_type(get<0>(p1))) + // * (return_type(get<1>(p2)) + return_type(get<1>(p1))); + + // SUM += (x1 + x2) * (y1 - y2) + state.sum += (return_type(get<0>(p1)) + return_type(get<0>(p2))) + * (return_type(get<1>(p1)) - return_type(get<1>(p2))); } static inline return_type result(summation const& state) diff --git a/boost/geometry/strategies/cartesian/box_in_box.hpp b/boost/geometry/strategies/cartesian/box_in_box.hpp index 56aef9e4df..fe77de9e2d 100644 --- a/boost/geometry/strategies/cartesian/box_in_box.hpp +++ b/boost/geometry/strategies/cartesian/box_in_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 2015. -// Modifications copyright (c) 2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2015, 2016. +// Modifications copyright (c) 2016, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -25,6 +25,7 @@ #include <boost/geometry/core/coordinate_dimension.hpp> #include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/within.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> namespace boost { namespace geometry { namespace strategy @@ -34,6 +35,8 @@ namespace boost { namespace geometry { namespace strategy namespace within { + +template <typename Geometry, std::size_t Dimension, typename CSTag> struct box_within_range { template <typename BoxContainedValue, typename BoxContainingValue> @@ -48,6 +51,7 @@ struct box_within_range }; +template <typename Geometry, std::size_t Dimension, typename CSTag> struct box_covered_by_range { template <typename BoxContainedValue, typename BoxContainingValue> @@ -61,9 +65,84 @@ struct box_covered_by_range }; +struct box_within_longitude_check +{ + template <typename CalcT> + static inline bool apply(CalcT const& diff_ed) + { + return diff_ed > CalcT(0); + } +}; + +struct box_covered_by_longitude_check +{ + template <typename CalcT> + static inline bool apply(CalcT const&) + { + return true; + } +}; + +template <typename Geometry, + typename InteriorCheck> +struct box_longitude_range +{ + template <typename BoxContainedValue, typename BoxContainingValue> + static inline bool apply(BoxContainedValue const& bed_min, + BoxContainedValue const& bed_max, + BoxContainingValue const& bing_min, + BoxContainingValue const& bing_max) + { + typedef typename select_most_precise + < + BoxContainedValue, + BoxContainingValue + >::type calc_t; + typedef typename coordinate_system<Geometry>::type::units units_t; + typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; + + // min <= max <=> diff >= 0 + calc_t const diff_ed = bed_max - bed_min; + calc_t const diff_ing = bing_max - bing_min; + + // if containing covers the whole globe it contains all + if (diff_ing >= constants::period()) + { + return true; + } + + // if containing is smaller it cannot contain + // and check interior (within vs covered_by) + if (diff_ing < diff_ed || ! InteriorCheck::apply(diff_ed)) + { + return false; + } + + // calculate positive longitude translation with bing_min as origin + calc_t const diff_min = math::longitude_distance_unsigned<units_t>(bing_min, bed_min); + + // max of contained translated into the containing origin must be lesser than max of containing + return bing_min + diff_min + diff_ed <= bing_max; + } +}; + + +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename Geometry> +struct box_within_range<Geometry, 0, spherical_tag> + : box_longitude_range<Geometry, box_within_longitude_check> +{}; + + +template <typename Geometry> +struct box_covered_by_range<Geometry, 0, spherical_tag> + : box_longitude_range<Geometry, box_covered_by_longitude_check> +{}; + + template < - typename SubStrategy, + template <typename, std::size_t, typename> class SubStrategy, typename Box1, typename Box2, std::size_t Dimension, @@ -74,8 +153,9 @@ struct relate_box_box_loop static inline bool apply(Box1 const& b_contained, Box2 const& b_containing) { assert_dimension_equal<Box1, Box2>(); + typedef typename tag_cast<typename cs_tag<Box1>::type, spherical_tag>::type cs_tag_t; - if (! SubStrategy::apply( + if (! SubStrategy<Box1, Dimension, cs_tag_t>::apply( get<min_corner, Dimension>(b_contained), get<max_corner, Dimension>(b_contained), get<min_corner, Dimension>(b_containing), @@ -97,7 +177,7 @@ struct relate_box_box_loop template < - typename SubStrategy, + template <typename, std::size_t, typename> class SubStrategy, typename Box1, typename Box2, std::size_t DimensionCount @@ -114,7 +194,7 @@ template < typename Box1, typename Box2, - typename SubStrategy = box_within_range + template <typename, std::size_t, typename> class SubStrategy = box_within_range > struct box_in_box { @@ -150,6 +230,19 @@ struct default_strategy typedef within::box_in_box<BoxContained, BoxContaining> type; }; +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename BoxContained, typename BoxContaining> +struct default_strategy + < + box_tag, box_tag, + box_tag, areal_tag, + spherical_tag, spherical_tag, + BoxContained, BoxContaining + > +{ + typedef within::box_in_box<BoxContained, BoxContaining> type; +}; + }} // namespace within::services @@ -172,6 +265,24 @@ struct default_strategy > type; }; +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename BoxContained, typename BoxContaining> +struct default_strategy + < + box_tag, box_tag, + box_tag, areal_tag, + spherical_tag, spherical_tag, + BoxContained, BoxContaining + > +{ + typedef within::box_in_box + < + BoxContained, BoxContaining, + within::box_covered_by_range + > type; +}; + + }} // namespace covered_by::services diff --git a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp index 47763a9f30..d081173a3e 100644 --- a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp +++ b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @@ -25,7 +25,7 @@ #include <boost/math/special_functions/fpclassify.hpp> #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/cast.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/core/coordinate_type.hpp> diff --git a/boost/geometry/strategies/cartesian/distance_projected_point.hpp b/boost/geometry/strategies/cartesian/distance_projected_point.hpp index dff4a77f6f..b7127c834d 100644 --- a/boost/geometry/strategies/cartesian/distance_projected_point.hpp +++ b/boost/geometry/strategies/cartesian/distance_projected_point.hpp @@ -22,7 +22,7 @@ #include <boost/concept_check.hpp> #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/point_type.hpp> diff --git a/boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp b/boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp index c4090044f4..7ff86c7cfa 100644 --- a/boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp +++ b/boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp @@ -24,8 +24,6 @@ #include <algorithm> #include <boost/concept_check.hpp> -#include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/point_type.hpp> diff --git a/boost/geometry/strategies/cartesian/distance_pythagoras.hpp b/boost/geometry/strategies/cartesian/distance_pythagoras.hpp index 77bdc966a6..665426ecb0 100644 --- a/boost/geometry/strategies/cartesian/distance_pythagoras.hpp +++ b/boost/geometry/strategies/cartesian/distance_pythagoras.hpp @@ -15,9 +15,6 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_DISTANCE_PYTHAGORAS_HPP -#include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> - #include <boost/geometry/core/access.hpp> #include <boost/geometry/strategies/distance.hpp> @@ -26,8 +23,6 @@ #include <boost/geometry/util/calculation_type.hpp> - - namespace boost { namespace geometry { diff --git a/boost/geometry/strategies/cartesian/point_in_box.hpp b/boost/geometry/strategies/cartesian/point_in_box.hpp index bd2303cbc4..28ed8215e7 100644 --- a/boost/geometry/strategies/cartesian/point_in_box.hpp +++ b/boost/geometry/strategies/cartesian/point_in_box.hpp @@ -4,8 +4,8 @@ // 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 Adam Wulkiewicz, on behalf of Oracle @@ -24,6 +24,7 @@ #include <boost/geometry/core/coordinate_dimension.hpp> #include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/within.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> namespace boost { namespace geometry { namespace strategy @@ -33,6 +34,7 @@ namespace within { +template <typename Geometry, std::size_t Dimension, typename CSTag> struct within_range { template <typename Value1, typename Value2> @@ -43,6 +45,7 @@ struct within_range }; +template <typename Geometry, std::size_t Dimension, typename CSTag> struct covered_by_range { template <typename Value1, typename Value2> @@ -53,9 +56,76 @@ struct covered_by_range }; +// NOTE: the result would be the same if instead of structs defined below +// the above xxx_range were used with the following arguments: +// (min_value + diff_min, min_value, max_value) +struct within_longitude_range +{ + template <typename CalcT> + static inline bool apply(CalcT const& diff_min, CalcT const& min_value, CalcT const& max_value) + { + CalcT const c0 = 0; + return diff_min > c0 && min_value + diff_min < max_value; + } +}; + +struct covered_by_longitude_range +{ + template <typename CalcT> + static inline bool apply(CalcT const& diff_min, CalcT const& min_value, CalcT const& max_value) + { + return min_value + diff_min <= max_value; + } +}; + + +template <typename Geometry, + typename ResultCheck> +struct longitude_range +{ + template <typename Value1, typename Value2> + static inline bool apply(Value1 const& value, Value2 const& min_value, Value2 const& max_value) + { + typedef typename select_most_precise + < + Value1, Value2 + >::type calc_t; + typedef typename coordinate_system<Geometry>::type::units units_t; + typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; + + // min <= max <=> diff >= 0 + calc_t const diff_ing = max_value - min_value; + + // if containing covers the whole globe it contains all + if (diff_ing >= constants::period()) + { + return true; + } + + // calculate positive longitude translation with min_value as origin + calc_t const diff_min = math::longitude_distance_unsigned<units_t, calc_t>(min_value, value); + + return ResultCheck::template apply<calc_t>(diff_min, min_value, max_value); + } +}; + + +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename Geometry> +struct within_range<Geometry, 0, spherical_tag> + : longitude_range<Geometry, within_longitude_range> +{}; + + +template <typename Geometry> +struct covered_by_range<Geometry, 0, spherical_tag> + : longitude_range<Geometry, covered_by_longitude_range> +{}; + + template < - typename SubStrategy, + template <typename, std::size_t, typename> class SubStrategy, typename Point, typename Box, std::size_t Dimension, @@ -65,7 +135,9 @@ struct relate_point_box_loop { static inline bool apply(Point const& point, Box const& box) { - if (! SubStrategy::apply(get<Dimension>(point), + typedef typename tag_cast<typename cs_tag<Point>::type, spherical_tag>::type cs_tag_t; + + if (! SubStrategy<Point, Dimension, cs_tag_t>::apply(get<Dimension>(point), get<min_corner, Dimension>(box), get<max_corner, Dimension>(box)) ) @@ -85,7 +157,7 @@ struct relate_point_box_loop template < - typename SubStrategy, + template <typename, std::size_t, typename> class SubStrategy, typename Point, typename Box, std::size_t DimensionCount @@ -103,7 +175,7 @@ template < typename Point, typename Box, - typename SubStrategy = within_range + template <typename, std::size_t, typename> class SubStrategy = within_range > struct point_in_box { @@ -140,6 +212,19 @@ struct default_strategy typedef within::point_in_box<Point, Box> type; }; +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename Point, typename Box> +struct default_strategy + < + point_tag, box_tag, + point_tag, areal_tag, + spherical_tag, spherical_tag, + Point, Box + > +{ + typedef within::point_in_box<Point, Box> type; +}; + }} // namespace within::services @@ -164,6 +249,23 @@ struct default_strategy > type; }; +// spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag +template <typename Point, typename Box> +struct default_strategy + < + point_tag, box_tag, + point_tag, areal_tag, + spherical_tag, spherical_tag, + Point, Box + > +{ + typedef within::point_in_box + < + Point, Box, + within::covered_by_range + > type; +}; + }} // namespace covered_by::services diff --git a/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/boost/geometry/strategies/cartesian/side_by_triangle.hpp index 77443d46a9..ba7749ba76 100644 --- a/boost/geometry/strategies/cartesian/side_by_triangle.hpp +++ b/boost/geometry/strategies/cartesian/side_by_triangle.hpp @@ -20,7 +20,8 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_BY_TRIANGLE_HPP #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/core/access.hpp> diff --git a/boost/geometry/strategies/geographic/distance_andoyer.hpp b/boost/geometry/strategies/geographic/distance_andoyer.hpp index 64de8c1a41..1646727d09 100644 --- a/boost/geometry/strategies/geographic/distance_andoyer.hpp +++ b/boost/geometry/strategies/geographic/distance_andoyer.hpp @@ -1,9 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-2016 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016. +// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -20,6 +20,7 @@ #include <boost/geometry/core/radius.hpp> #include <boost/geometry/core/srs.hpp> +#include <boost/geometry/algorithms/detail/andoyer_inverse.hpp> #include <boost/geometry/algorithms/detail/flattening.hpp> #include <boost/geometry/strategies/distance.hpp> @@ -50,6 +51,8 @@ are about the same as Vincenty. In my (Barend's) testcases the results didn't di \see http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115 (implementation) \see http://futureboy.homeip.net/frinksamp/navigation.frink (implementation) \see http://www.voidware.com/earthdist.htm (implementation) +\see http://www.dtic.mil/docs/citations/AD0627893 +\see http://www.dtic.mil/docs/citations/AD703541 */ template < @@ -82,16 +85,17 @@ public : : m_spheroid(spheroid) {} - template <typename Point1, typename Point2> inline typename calculation_type<Point1, Point2>::type apply(Point1 const& point1, Point2 const& point2) const { - return calc<typename calculation_type<Point1, Point2>::type> - ( - get_as_radian<0>(point1), get_as_radian<1>(point1), - get_as_radian<0>(point2), get_as_radian<1>(point2) - ); + return geometry::detail::andoyer_inverse + < + typename calculation_type<Point1, Point2>::type, + true, false + >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1), + get_as_radian<0>(point2), get_as_radian<1>(point2), + m_spheroid).distance; } inline Spheroid const& model() const @@ -100,55 +104,6 @@ public : } private : - template <typename CT, typename T> - inline CT calc(T const& lon1, - T const& lat1, - T const& lon2, - T const& lat2) const - { - CT const G = (lat1 - lat2) / 2.0; - CT const lambda = (lon1 - lon2) / 2.0; - - if (geometry::math::equals(lambda, 0.0) - && geometry::math::equals(G, 0.0)) - { - return 0.0; - } - - CT const F = (lat1 + lat2) / 2.0; - - CT const sinG2 = math::sqr(sin(G)); - CT const cosG2 = math::sqr(cos(G)); - CT const sinF2 = math::sqr(sin(F)); - CT const cosF2 = math::sqr(cos(F)); - CT const sinL2 = math::sqr(sin(lambda)); - CT const cosL2 = math::sqr(cos(lambda)); - - CT const S = sinG2 * cosL2 + cosF2 * sinL2; - CT const C = cosG2 * cosL2 + sinF2 * sinL2; - - CT const c0 = 0; - CT const c1 = 1; - CT const c2 = 2; - CT const c3 = 3; - - if (geometry::math::equals(S, c0) || geometry::math::equals(C, c0)) - { - return c0; - } - - CT const radius_a = CT(get_radius<0>(m_spheroid)); - CT const flattening = geometry::detail::flattening<CT>(m_spheroid); - - CT const omega = atan(math::sqrt(S / C)); - CT const r3 = c3 * math::sqrt(S * C) / omega; // not sure if this is r or greek nu - CT const D = c2 * omega * radius_a; - CT const H1 = (r3 - c1) / (c2 * C); - CT const H2 = (r3 + c1) / (c2 * S); - - return D * (c1 + flattening * (H1 * sinF2 * cosG2 - H2 * cosF2 * sinG2) ); - } - Spheroid m_spheroid; }; diff --git a/boost/geometry/strategies/spherical/ssf.hpp b/boost/geometry/strategies/spherical/ssf.hpp index 81f3205e90..c99e5835e2 100644 --- a/boost/geometry/strategies/spherical/ssf.hpp +++ b/boost/geometry/strategies/spherical/ssf.hpp @@ -2,6 +2,10 @@ // Copyright (c) 2011-2012 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) @@ -9,8 +13,6 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP -#include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/access.hpp> @@ -61,9 +63,9 @@ int spherical_side_formula(T const& lambda1, T const& delta1, + (c1x * c2y - c1y * c2x) * sin(delta); T zero = T(); - return dist > zero ? 1 - : dist < zero ? -1 - : 0; + return math::equals(dist, zero) ? 0 + : dist > zero ? 1 + : -1; // dist < zero } } diff --git a/boost/geometry/util/calculation_type.hpp b/boost/geometry/util/calculation_type.hpp index 2ac0de87d9..18eac4fbb7 100644 --- a/boost/geometry/util/calculation_type.hpp +++ b/boost/geometry/util/calculation_type.hpp @@ -13,7 +13,9 @@ #include <boost/config.hpp> #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/type_traits/is_fundamental.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/util/select_coordinate_type.hpp> #include <boost/geometry/util/select_most_precise.hpp> diff --git a/boost/geometry/util/combine_if.hpp b/boost/geometry/util/combine_if.hpp index cfc3449267..5d94c34461 100644 --- a/boost/geometry/util/combine_if.hpp +++ b/boost/geometry/util/combine_if.hpp @@ -24,8 +24,6 @@ #include <boost/mpl/insert.hpp> #include <boost/mpl/placeholders.hpp> -#include <boost/type_traits.hpp> - namespace boost { namespace geometry { diff --git a/boost/geometry/util/math.hpp b/boost/geometry/util/math.hpp index c193c8f3f1..234cfa1ed0 100644 --- a/boost/geometry/util/math.hpp +++ b/boost/geometry/util/math.hpp @@ -32,6 +32,8 @@ #include <boost/type_traits/is_fundamental.hpp> #include <boost/type_traits/is_integral.hpp> +#include <boost/geometry/core/cs.hpp> + #include <boost/geometry/util/select_most_precise.hpp> namespace boost { namespace geometry @@ -601,6 +603,65 @@ inline T r2d() } +#ifndef DOXYGEN_NO_DETAIL +namespace detail { + +template <typename DegreeOrRadian> +struct as_radian +{ + template <typename T> + static inline T apply(T const& value) + { + return value; + } +}; + +template <> +struct as_radian<degree> +{ + template <typename T> + static inline T apply(T const& value) + { + return value * d2r<T>(); + } +}; + +template <typename DegreeOrRadian> +struct from_radian +{ + template <typename T> + static inline T apply(T const& value) + { + return value; + } +}; + +template <> +struct from_radian<degree> +{ + template <typename T> + static inline T apply(T const& value) + { + return value * r2d<T>(); + } +}; + +} // namespace detail +#endif + +template <typename DegreeOrRadian, typename T> +inline T as_radian(T const& value) +{ + return detail::as_radian<DegreeOrRadian>::apply(value); +} + +template <typename DegreeOrRadian, typename T> +inline T from_radian(T const& value) +{ + return detail::from_radian<DegreeOrRadian>::apply(value); +} + + /*! \brief Calculates the haversine of an angle \ingroup utility diff --git a/boost/geometry/util/normalize_spheroidal_coordinates.hpp b/boost/geometry/util/normalize_spheroidal_coordinates.hpp index 72dfa5a6f1..377051eef1 100644 --- a/boost/geometry/util/normalize_spheroidal_coordinates.hpp +++ b/boost/geometry/util/normalize_spheroidal_coordinates.hpp @@ -1,8 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2015, Oracle and/or its affiliates. +// 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 // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -119,6 +120,27 @@ protected: } public: + static inline void apply(CoordinateType& longitude) + { + // normalize longitude + if (math::equals(math::abs(longitude), constants::half_period())) + { + longitude = constants::half_period(); + } + else if (longitude > constants::half_period()) + { + longitude = normalize_up(longitude); + if (math::equals(longitude, -constants::half_period())) + { + longitude = constants::half_period(); + } + } + else if (longitude < -constants::half_period()) + { + longitude = normalize_down(longitude); + } + } + static inline void apply(CoordinateType& longitude, CoordinateType& latitude, bool normalize_poles = true) @@ -148,22 +170,7 @@ public: #endif // BOOST_GEOMETRY_NORMALIZE_LATITUDE // normalize longitude - if (math::equals(math::abs(longitude), constants::half_period())) - { - longitude = constants::half_period(); - } - else if (longitude > constants::half_period()) - { - longitude = normalize_up(longitude); - if (math::equals(longitude, -constants::half_period())) - { - longitude = constants::half_period(); - } - } - else if (longitude < -constants::half_period()) - { - longitude = normalize_down(longitude); - } + apply(longitude); // finally normalize poles if (normalize_poles) @@ -210,6 +217,71 @@ inline void normalize_spheroidal_coordinates(CoordinateType& longitude, } +/*! +\brief Short utility to normalize the longitude on a spheroid. + Note that in general both coordinates should be normalized at once. + This utility is suitable e.g. for normalization of the difference of longitudes. +\tparam Units The units of the coordindate system in the spheroid +\tparam CoordinateType The type of the coordinates +\param longitude Longitude +\ingroup utility +*/ +template <typename Units, typename CoordinateType> +inline void normalize_longitude(CoordinateType& longitude) +{ + detail::normalize_spheroidal_coordinates + < + Units, CoordinateType + >::apply(longitude); +} + + +/*! +\brief Short utility to calculate difference between two longitudes + normalized in range (-180, 180]. +\tparam Units The units of the coordindate system in the spheroid +\tparam CoordinateType The type of the coordinates +\param longitude1 Longitude 1 +\param longitude2 Longitude 2 +\ingroup utility +*/ +template <typename Units, typename CoordinateType> +inline CoordinateType longitude_distance_signed(CoordinateType const& longitude1, + CoordinateType const& longitude2) +{ + CoordinateType diff = longitude2 - longitude1; + math::normalize_longitude<Units, CoordinateType>(diff); + return diff; +} + + +/*! +\brief Short utility to calculate difference between two longitudes + normalized in range [0, 360). +\tparam Units The units of the coordindate system in the spheroid +\tparam CoordinateType The type of the coordinates +\param longitude1 Longitude 1 +\param longitude2 Longitude 2 +\ingroup utility +*/ +template <typename Units, typename CoordinateType> +inline CoordinateType longitude_distance_unsigned(CoordinateType const& longitude1, + CoordinateType const& longitude2) +{ + typedef math::detail::constants_on_spheroid + < + CoordinateType, Units + > constants; + + CoordinateType const c0 = 0; + CoordinateType diff = longitude_distance_signed<Units>(longitude1, longitude2); + if (diff < c0) // (-180, 180] -> [0, 360) + { + diff += constants::period(); + } + return diff; +} + } // namespace math diff --git a/boost/geometry/util/parameter_type_of.hpp b/boost/geometry/util/parameter_type_of.hpp index 3ed09ab9bb..731fb55a56 100644 --- a/boost/geometry/util/parameter_type_of.hpp +++ b/boost/geometry/util/parameter_type_of.hpp @@ -21,7 +21,7 @@ #include <boost/mpl/at.hpp> #include <boost/mpl/int.hpp> #include <boost/mpl/plus.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/remove_reference.hpp> namespace boost { namespace geometry diff --git a/boost/geometry/util/promote_floating_point.hpp b/boost/geometry/util/promote_floating_point.hpp index 0c74cb8d6c..b84069db1b 100644 --- a/boost/geometry/util/promote_floating_point.hpp +++ b/boost/geometry/util/promote_floating_point.hpp @@ -16,7 +16,7 @@ #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_integral.hpp> namespace boost { namespace geometry diff --git a/boost/geometry/util/range.hpp b/boost/geometry/util/range.hpp index 0f480df58a..b10a008059 100644 --- a/boost/geometry/util/range.hpp +++ b/boost/geometry/util/range.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 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, 2014, 2015, 2016. +// Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -15,14 +15,21 @@ #define BOOST_GEOMETRY_UTIL_RANGE_HPP #include <algorithm> +#include <iterator> #include <boost/concept_check.hpp> #include <boost/config.hpp> +#include <boost/core/addressof.hpp> #include <boost/range/concepts.hpp> #include <boost/range/begin.hpp> #include <boost/range/end.hpp> #include <boost/range/empty.hpp> +#include <boost/range/difference_type.hpp> +#include <boost/range/iterator.hpp> +#include <boost/range/rbegin.hpp> +#include <boost/range/reference.hpp> #include <boost/range/size.hpp> +#include <boost/range/value_type.hpp> #include <boost/type_traits/is_convertible.hpp> #include <boost/geometry/core/assert.hpp> @@ -362,6 +369,50 @@ erase(Range & rng, return erase(rng, first, last); } +// back_inserter + +template <class Container> +class back_insert_iterator + : public std::iterator<std::output_iterator_tag, void, void, void, void> +{ +public: + typedef Container container_type; + + explicit back_insert_iterator(Container & c) + : container(boost::addressof(c)) + {} + + back_insert_iterator & operator=(typename Container::value_type const& value) + { + range::push_back(*container, value); + return *this; + } + + back_insert_iterator & operator* () + { + return *this; + } + + back_insert_iterator & operator++ () + { + return *this; + } + + back_insert_iterator operator++(int) + { + return *this; + } + +private: + Container * container; +}; + +template <typename Range> +inline back_insert_iterator<Range> back_inserter(Range & rng) +{ + return back_insert_iterator<Range>(rng); +} + }}} // namespace boost::geometry::range #endif // BOOST_GEOMETRY_UTIL_RANGE_HPP diff --git a/boost/geometry/util/select_calculation_type.hpp b/boost/geometry/util/select_calculation_type.hpp index 0e3faf229b..b6405831eb 100644 --- a/boost/geometry/util/select_calculation_type.hpp +++ b/boost/geometry/util/select_calculation_type.hpp @@ -21,7 +21,7 @@ #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/util/select_coordinate_type.hpp> diff --git a/boost/geometry/util/select_most_precise.hpp b/boost/geometry/util/select_most_precise.hpp index b5b1388cf8..3169f5fa17 100644 --- a/boost/geometry/util/select_most_precise.hpp +++ b/boost/geometry/util/select_most_precise.hpp @@ -20,7 +20,8 @@ #define BOOST_GEOMETRY_UTIL_SELECT_MOST_PRECISE_HPP #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/type_traits/is_fundamental.hpp> namespace boost { namespace geometry |