diff options
Diffstat (limited to 'boost/geometry/algorithms/detail/buffer')
8 files changed, 1357 insertions, 393 deletions
diff --git a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp index c959ee849b..127e4c3fb2 100644 --- a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp @@ -12,6 +12,7 @@ #include <cstddef> #include <iterator> +#include <boost/core/ignore_unused.hpp> #include <boost/numeric/conversion/cast.hpp> #include <boost/range.hpp> @@ -20,6 +21,7 @@ #include <boost/geometry/core/exterior_ring.hpp> #include <boost/geometry/core/interior_rings.hpp> +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/strategies/buffer.hpp> @@ -28,6 +30,7 @@ #include <boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp> #include <boost/geometry/algorithms/detail/buffer/parallel_continue.hpp> +#include <boost/geometry/algorithms/num_interior_rings.hpp> #include <boost/geometry/algorithms/simplify.hpp> #include <boost/geometry/views/detail/normalized_view.hpp> @@ -144,15 +147,25 @@ struct buffer_range intersection_point); } + switch(join) { case strategy::buffer::join_continue : // No join, we get two consecutive sides - return; + break; case strategy::buffer::join_concave : - collection.add_piece(strategy::buffer::buffered_concave, - previous_input, prev_perp2, perp1); - return; + { + std::vector<output_point_type> range_out; + range_out.push_back(prev_perp2); + range_out.push_back(previous_input); + collection.add_piece(strategy::buffer::buffered_concave, previous_input, range_out); + + range_out.clear(); + range_out.push_back(previous_input); + range_out.push_back(perp1); + collection.add_piece(strategy::buffer::buffered_concave, previous_input, range_out); + } + break; case strategy::buffer::join_spike : { // For linestrings, only add spike at one side to avoid @@ -160,22 +173,24 @@ struct buffer_range std::vector<output_point_type> range_out; end_strategy.apply(penultimate_input, prev_perp2, previous_input, perp1, side, distance, range_out); collection.add_endcap(end_strategy, range_out, previous_input); + collection.set_current_ring_concave(); } - return; + break; case strategy::buffer::join_convex : - break; // All code below handles this - } - - // The corner is convex, we create a join - // TODO (future) - avoid a separate vector, add the piece directly - std::vector<output_point_type> range_out; - if (join_strategy.apply(intersection_point, - previous_input, prev_perp2, perp1, - distance.apply(previous_input, input, side), - range_out)) - { - collection.add_piece(strategy::buffer::buffered_join, - previous_input, range_out); + { + // The corner is convex, we create a join + // TODO (future) - avoid a separate vector, add the piece directly + std::vector<output_point_type> range_out; + if (join_strategy.apply(intersection_point, + previous_input, prev_perp2, perp1, + distance.apply(previous_input, input, side), + range_out)) + { + collection.add_piece(strategy::buffer::buffered_join, + previous_input, range_out); + } + } + break; } } @@ -225,18 +240,13 @@ struct buffer_range output_point_type& last_p1, output_point_type& last_p2) { + boost::ignore_unused(side_strategy); + typedef typename std::iterator_traits < Iterator >::value_type point_type; - typedef typename robust_point_type - < - point_type, - RobustPolicy - >::type robust_point_type; - - robust_point_type previous_robust_input; point_type second_point, penultimate_point, ultimate_point; // last two points from begin/end /* @@ -261,57 +271,50 @@ struct buffer_range Iterator it = begin; - geometry::recalculate(previous_robust_input, *begin, robust_policy); - std::vector<output_point_type> generated_side; generated_side.reserve(2); for (Iterator prev = it++; it != end; ++it) { - robust_point_type robust_input; - geometry::recalculate(robust_input, *it, robust_policy); - // Check on equality - however, if input is simplified, this is - // unlikely (though possible by rescaling or for degenerated pointlike polygons) - if (! detail::equals::equals_point_point(previous_robust_input, robust_input)) - { - generated_side.clear(); - side_strategy.apply(*prev, *it, side, - distance_strategy, generated_side); + generated_side.clear(); + side_strategy.apply(*prev, *it, side, + distance_strategy, generated_side); - if (generated_side.empty()) - { - break; - } + if (generated_side.empty()) + { + // Because input is simplified, this is improbable, + // but it can happen for degenerate geometries + // Further handling of this side is skipped + continue; + } - result = true; + result = true; - if (! first) - { - add_join(collection, - penultimate_point, - *prev, last_p1, last_p2, - *it, generated_side.front(), generated_side.back(), - side, - distance_strategy, join_strategy, end_strategy, - robust_policy); - } + if (! first) + { + add_join(collection, + penultimate_point, + *prev, last_p1, last_p2, + *it, generated_side.front(), generated_side.back(), + side, + distance_strategy, join_strategy, end_strategy, + robust_policy); + } - collection.add_side_piece(*prev, *it, generated_side, first); + collection.add_side_piece(*prev, *it, generated_side, first); - penultimate_point = *prev; - ultimate_point = *it; - last_p1 = generated_side.front(); - last_p2 = generated_side.back(); - prev = it; - if (first) - { - first = false; - second_point = *it; - first_p1 = generated_side.front(); - first_p2 = generated_side.back(); - } + penultimate_point = *prev; + ultimate_point = *it; + last_p1 = generated_side.front(); + last_p2 = generated_side.back(); + prev = it; + if (first) + { + first = false; + second_point = *it; + first_p1 = generated_side.front(); + first_p2 = generated_side.back(); } - previous_robust_input = robust_input; } return result; } @@ -775,7 +778,7 @@ public: distance, side_strategy, join_strategy, end_strategy, point_strategy, robust_policy); - collection.finish_ring(); + collection.finish_ring(false, geometry::num_interior_rings(polygon) > 0u); } apply_interior_rings(interior_rings(polygon), @@ -839,6 +842,8 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator VisitPiecesPolicy& visit_pieces_policy ) { + boost::ignore_unused(visit_pieces_policy); + typedef detail::buffer::buffered_piece_collection < typename geometry::ring_type<GeometryOutput>::type, @@ -852,6 +857,11 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator typename tag_cast<typename tag<GeometryInput>::type, areal_tag>::type, areal_tag >::type::value; + bool const linear = boost::is_same + < + typename tag_cast<typename tag<GeometryInput>::type, linear_tag>::type, + linear_tag + >::type::value; dispatch::buffer_inserter < @@ -868,9 +878,10 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator robust_policy); collection.get_turns(); - if (areal) + collection.classify_turns(linear); + if (BOOST_GEOMETRY_CONDITION(areal)) { - collection.check_remaining_points(distance_strategy.factor()); + collection.check_remaining_points(distance_strategy); } // Visit the piece collection. This does nothing (by default), but @@ -889,7 +900,8 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator // - the output is counter clockwise // and avoid reversing twice bool reverse = distance_strategy.negative() && areal; - if (geometry::point_order<GeometryOutput>::value == counterclockwise) + if (BOOST_GEOMETRY_CONDITION( + geometry::point_order<GeometryOutput>::value == counterclockwise)) { reverse = ! reverse; } @@ -898,6 +910,11 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator collection.reverse(); } + if (distance_strategy.negative() && areal) + { + collection.discard_nonintersecting_deflated_rings(); + } + collection.template assign<GeometryOutput>(out); // Visit collection again diff --git a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp index 6a2e6b32c5..c5bb8acc05 100644 --- a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp @@ -35,7 +35,7 @@ namespace detail { namespace buffer enum intersection_location_type { - location_ok, inside_buffer, inside_original + location_ok, inside_buffer, location_discard }; class backtrack_for_buffer @@ -120,10 +120,14 @@ struct buffer_turn_info return robust_point; } - intersection_location_type location; int count_within; + + bool within_original; + int count_on_original_boundary; + int count_in_original; // increased by +1 for in ext.ring, -1 for int.ring + int count_on_offsetted; int count_on_helper; int count_within_near_offsetted; @@ -138,6 +142,9 @@ struct buffer_turn_info : turn_index(-1) , location(location_ok) , count_within(0) + , within_original(false) + , count_on_original_boundary(0) + , count_in_original(0) , count_on_offsetted(0) , count_on_helper(0) , count_within_near_offsetted(0) diff --git a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp index 558a61fcb4..a501e3f197 100644 --- a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp @@ -12,8 +12,9 @@ #include <algorithm> #include <cstddef> #include <set> -#include <boost/range.hpp> +#include <boost/core/ignore_unused.hpp> +#include <boost/range.hpp> #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/point_type.hpp> @@ -24,13 +25,14 @@ #include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/geometries/ring.hpp> -#include <boost/geometry/geometries/polygon.hpp> #include <boost/geometry/algorithms/detail/buffer/buffered_ring.hpp> #include <boost/geometry/algorithms/detail/buffer/buffer_policies.hpp> #include <boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp> #include <boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp> +#include <boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp> +#include <boost/geometry/algorithms/detail/disjoint/point_box.hpp> #include <boost/geometry/algorithms/detail/overlay/add_rings.hpp> #include <boost/geometry/algorithms/detail/overlay/assign_parents.hpp> #include <boost/geometry/algorithms/detail/overlay/enrichment_info.hpp> @@ -41,6 +43,8 @@ #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> #include <boost/geometry/algorithms/detail/occupation_info.hpp> #include <boost/geometry/algorithms/detail/partition.hpp> +#include <boost/geometry/algorithms/detail/sections/sectionalize.hpp> +#include <boost/geometry/algorithms/detail/sections/section_box_policies.hpp> #include <boost/geometry/util/range.hpp> @@ -89,7 +93,7 @@ enum segment_relation_code * form together the offsetted ring (marked with o below) * The 8 pieces are part of the piece collection and use for inside-checks * The inner parts form (using 1 or 2 points per piece, often co-located) - * form together the robust_ring (marked with r below) + * form together the robust_polygons (marked with r below) * The remaining piece-segments are helper-segments (marked with h) * * ooooooooooooooooo @@ -120,7 +124,7 @@ struct buffered_piece_collection // Robust ring/polygon type, always clockwise typedef geometry::model::ring<robust_point_type> robust_ring_type; - typedef geometry::model::polygon<robust_point_type> robust_polygon_type; + typedef geometry::model::box<robust_point_type> robust_box_type; typedef typename strategy::side::services::default_strategy < @@ -164,6 +168,9 @@ struct buffered_piece_collection struct piece { + typedef robust_ring_type piece_robust_ring_type; + typedef geometry::section<robust_box_type, 1> section_type; + strategy::buffer::piece_type type; int index; @@ -181,18 +188,58 @@ struct buffered_piece_collection #if defined(BOOST_GEOMETRY_BUFFER_USE_HELPER_POINTS) // 2: half, not part of offsetted rings - part of robust ring - std::vector<point_type> helper_points; // 4 points for segment, 3 points for join - 0 points for flat-end + std::vector<point_type> helper_points; // 4 points for side, 3 points for join - 0 points for flat-end #endif + bool is_monotonic_increasing[2]; // 0=x, 1=y + bool is_monotonic_decreasing[2]; // 0=x, 1=y + + // Monotonic sections of pieces around points + std::vector<section_type> sections; + + // Robust representations // 3: complete ring robust_ring_type robust_ring; - geometry::model::box<robust_point_type> robust_envelope; + robust_box_type robust_envelope; + robust_box_type robust_offsetted_envelope; std::vector<robust_turn> robust_turns; // Used only in insert_rescaled_piece_turns - we might use a map instead }; + struct robust_original + { + typedef robust_ring_type original_robust_ring_type; + typedef geometry::sections<robust_box_type, 1> sections_type; + + inline robust_original() + : m_is_interior(false) + , m_has_interiors(true) + {} + + inline robust_original(robust_ring_type const& ring, + bool is_interior, bool has_interiors) + : m_ring(ring) + , m_is_interior(is_interior) + , m_has_interiors(has_interiors) + { + geometry::envelope(m_ring, m_box); + + // create monotonic sections in y-dimension + typedef boost::mpl::vector_c<std::size_t, 1> dimensions; + geometry::sectionalize<false, dimensions>(m_ring, + detail::no_rescale_policy(), m_sections); + } + + robust_ring_type m_ring; + robust_box_type m_box; + sections_type m_sections; + + bool m_is_interior; + bool m_has_interiors; + }; + typedef std::vector<piece> piece_vector_type; piece_vector_type m_pieces; @@ -200,11 +247,17 @@ struct buffered_piece_collection int m_first_piece_index; buffered_ring_collection<buffered_ring<Ring> > offsetted_rings; // indexed by multi_index - buffered_ring_collection<robust_polygon_type> robust_polygons; // robust representation of the original(s) + std::vector<robust_original> robust_originals; // robust representation of the original(s) robust_ring_type current_robust_ring; buffered_ring_collection<Ring> traversed_rings; segment_identifier current_segment_id; + // Specificly for offsetted rings around points + // but also for large joins with many points + typedef geometry::sections<robust_box_type, 2> sections_type; + sections_type monotonic_sections; + + RobustPolicy const& m_robust_policy; struct redundant_turn @@ -378,7 +431,7 @@ struct buffered_piece_collection } } - inline void classify_turns() + inline void classify_turns(bool linear) { for (typename boost::range_iterator<turn_vector_type>::type it = boost::begin(m_turns); it != boost::end(m_turns); ++it) @@ -387,6 +440,10 @@ struct buffered_piece_collection { it->location = inside_buffer; } + if (it->count_on_original_boundary > 0 && ! linear) + { + it->location = inside_buffer; + } if (it->count_within_near_offsetted > 0) { // Within can have in rare cases a rounding issue. We don't discard this @@ -394,31 +451,40 @@ struct buffered_piece_collection // will never start a new ring from this type of points. it->selectable_start = false; } - } } - inline void check_remaining_points(int factor) + template <typename DistanceStrategy> + inline void check_remaining_points(DistanceStrategy const& distance_strategy) { - // TODO: use partition + // Check if a turn is inside any of the originals + + turn_in_original_visitor<turn_vector_type> visitor(m_turns); + geometry::partition + < + robust_box_type, + turn_get_box, turn_in_original_ovelaps_box, + original_get_box, original_ovelaps_box, + include_turn_policy, detail::partition::include_all_policy + >::apply(m_turns, robust_originals, visitor); + + bool const deflate = distance_strategy.negative(); for (typename boost::range_iterator<turn_vector_type>::type it = boost::begin(m_turns); it != boost::end(m_turns); ++it) { - if (it->location == location_ok) + buffer_turn_info_type& turn = *it; + if (turn.location == location_ok) { - int code = -1; - for (std::size_t i = 0; i < robust_polygons.size(); i++) + if (deflate && turn.count_in_original <= 0) { - if (geometry::covered_by(it->robust_point, robust_polygons[i])) - { - code = 1; - break; - } + // For deflate: it is not in original, discard + turn.location = location_discard; } - if (code * factor == 1) + else if (! deflate && turn.count_in_original > 0) { - it->location = inside_original; + // For inflate: it is in original, discard + turn.location = location_discard; } } } @@ -474,6 +540,7 @@ struct buffered_piece_collection // Take into account for the box (intersection points should fall inside, // but in theory they can be one off because of rounding geometry::expand(pc.robust_envelope, it->robust_point); + geometry::expand(pc.robust_offsetted_envelope, it->robust_point); } } @@ -502,9 +569,10 @@ struct buffered_piece_collection { int const index_in_vector = 1 + rit->seg_id.segment_index - piece_segment_index; BOOST_ASSERT - ( - index_in_vector > 0 && index_in_vector < pc.offsetted_count - ); + ( + index_in_vector > 0 + && index_in_vector < pc.offsetted_count + ); pc.robust_ring.insert(boost::begin(pc.robust_ring) + index_in_vector, rit->point); pc.offsetted_count++; @@ -517,25 +585,132 @@ struct buffered_piece_collection BOOST_ASSERT(assert_indices_in_robust_rings()); } + template <std::size_t Dimension> + static inline void determine_monotonicity(piece& pc, + robust_point_type const& current, + robust_point_type const& next) + { + if (geometry::get<Dimension>(current) >= geometry::get<Dimension>(next)) + { + pc.is_monotonic_increasing[Dimension] = false; + } + if (geometry::get<Dimension>(current) <= geometry::get<Dimension>(next)) + { + pc.is_monotonic_decreasing[Dimension] = false; + } + } + + static inline void determine_properties(piece& pc) + { + pc.is_monotonic_increasing[0] = true; + pc.is_monotonic_increasing[1] = true; + pc.is_monotonic_decreasing[0] = true; + pc.is_monotonic_decreasing[1] = true; + + if (pc.offsetted_count < 2) + { + return; + } + + typename robust_ring_type::const_iterator current = pc.robust_ring.begin(); + typename robust_ring_type::const_iterator next = current + 1; + + for (int i = 1; i < pc.offsetted_count; i++) + { + determine_monotonicity<0>(pc, *current, *next); + determine_monotonicity<1>(pc, *current, *next); + current = next; + ++next; + } + + } + + void determine_properties() + { + for (typename piece_vector_type::iterator it = boost::begin(m_pieces); + it != boost::end(m_pieces); + ++it) + { + determine_properties(*it); + } + } + + inline void reverse_negative_robust_rings() + { + for (typename piece_vector_type::iterator it = boost::begin(m_pieces); + it != boost::end(m_pieces); + ++it) + { + piece& pc = *it; + if (geometry::area(pc.robust_ring) < 0) + { + // Rings can be ccw: + // - in a concave piece + // - in a line-buffer with a negative buffer-distance + std::reverse(pc.robust_ring.begin(), pc.robust_ring.end()); + } + } + } + + inline void prepare_buffered_point_piece(piece& pc) + { + // create monotonic sections in y-dimension + typedef boost::mpl::vector_c<std::size_t, 1> dimensions; + geometry::sectionalize<false, dimensions>(pc.robust_ring, + detail::no_rescale_policy(), pc.sections); + + // TODO (next phase) determine min/max radius + } + + inline void prepare_buffered_point_pieces() + { + for (typename piece_vector_type::iterator it = boost::begin(m_pieces); + it != boost::end(m_pieces); + ++it) + { + if (it->type == geometry::strategy::buffer::buffered_point) + { + prepare_buffered_point_piece(*it); + } + } + } + inline void get_turns() { + for(typename boost::range_iterator<sections_type>::type it + = boost::begin(monotonic_sections); + it != boost::end(monotonic_sections); + ++it) + { + enlarge_box(it->bounding_box, 1); + } + { // Calculate the turns piece_turn_visitor < + piece_vector_type, buffered_ring_collection<buffered_ring<Ring> >, turn_vector_type, RobustPolicy - > visitor(offsetted_rings, m_turns, m_robust_policy); + > visitor(m_pieces, offsetted_rings, m_turns, m_robust_policy); geometry::partition < - model::box<robust_point_type>, piece_get_box, piece_ovelaps_box - >::apply(m_pieces, visitor); + robust_box_type, + detail::section::get_section_box, + detail::section::overlaps_section_box + >::apply(monotonic_sections, visitor); } insert_rescaled_piece_turns(); + reverse_negative_robust_rings(); + + determine_properties(); + + prepare_buffered_point_pieces(); + { // Check if it is inside any of the pieces turn_in_piece_visitor @@ -545,16 +720,12 @@ struct buffered_piece_collection geometry::partition < - model::box<robust_point_type>, + robust_box_type, turn_get_box, turn_ovelaps_box, piece_get_box, piece_ovelaps_box >::apply(m_turns, m_pieces, visitor); } - - classify_turns(); - - //get_occupation(); } inline void start_new_ring() @@ -571,7 +742,38 @@ struct buffered_piece_collection m_first_piece_index = boost::size(m_pieces); } - inline void finish_ring(bool is_interior = false) + inline void update_closing_point() + { + BOOST_ASSERT(! offsetted_rings.empty()); + buffered_ring<Ring>& added = offsetted_rings.back(); + if (! boost::empty(added)) + { + range::back(added) = range::front(added); + } + } + + inline void update_last_point(point_type const& p, + buffered_ring<Ring>& ring) + { + // For the first point of a new piece, and there were already + // points in the offsetted ring, for some piece types the first point + // is a duplicate of the last point of the previous piece. + + // TODO: disable that, that point should not be added + + // For now, it is made equal because due to numerical instability, + // it can be a tiny bit off, possibly causing a self-intersection + + BOOST_ASSERT(boost::size(m_pieces) > 0); + if (! ring.empty() + && current_segment_id.segment_index + == m_pieces.back().first_seg_id.segment_index) + { + ring.back() = p; + } + } + + inline void finish_ring(bool is_interior = false, bool has_interiors = false) { if (m_first_piece_index == -1) { @@ -588,41 +790,50 @@ struct buffered_piece_collection } m_first_piece_index = -1; - if (!current_robust_ring.empty()) + update_closing_point(); + + if (! current_robust_ring.empty()) { - BOOST_ASSERT(geometry::equals(current_robust_ring.front(), current_robust_ring.back())); + BOOST_ASSERT + ( + geometry::equals(current_robust_ring.front(), + current_robust_ring.back()) + ); - if (is_interior) - { - if (!robust_polygons.empty()) - { - robust_polygons.back().inners().push_back(current_robust_ring); - } - } - else - { - robust_polygons.resize(robust_polygons.size() + 1); - robust_polygons.back().outer() = current_robust_ring; - } + robust_originals.push_back( + robust_original(current_robust_ring, + is_interior, has_interiors)); } } + inline void set_current_ring_concave() + { + BOOST_ASSERT(boost::size(offsetted_rings) > 0); + offsetted_rings.back().has_concave = true; + } + inline int add_point(point_type const& p) { - BOOST_ASSERT - ( - boost::size(offsetted_rings) > 0 - ); + BOOST_ASSERT(boost::size(offsetted_rings) > 0); + + buffered_ring<Ring>& current_ring = offsetted_rings.back(); + update_last_point(p, current_ring); current_segment_id.segment_index++; - offsetted_rings.back().push_back(p); - return offsetted_rings.back().size(); + current_ring.push_back(p); + return current_ring.size(); } //------------------------------------------------------------------------- - inline piece& create_piece(strategy::buffer::piece_type type, bool decrease_segment_index_by_one) + inline piece& create_piece(strategy::buffer::piece_type type, + bool decrease_segment_index_by_one) { + if (type == strategy::buffer::buffered_concave) + { + offsetted_rings.back().has_concave = true; + } + piece pc; pc.type = type; pc.index = boost::size(m_pieces); @@ -634,6 +845,7 @@ struct buffered_piece_collection std::size_t const n = boost::size(offsetted_rings.back()); pc.first_seg_id.segment_index = decrease_segment_index_by_one ? n - 1 : n; + pc.last_segment_index = pc.first_seg_id.segment_index; m_pieces.push_back(pc); return m_pieces.back(); @@ -641,6 +853,18 @@ struct buffered_piece_collection inline void init_rescale_piece(piece& pc, std::size_t helper_points_size) { + if (pc.first_seg_id.segment_index < 0) + { + // This indicates an error situation: an earlier piece was empty + // It currently does not happen + // std::cout << "EMPTY " << pc.type << " " << pc.index << " " << pc.first_seg_id.multi_index << std::endl; + pc.offsetted_count = 0; + return; + } + + BOOST_ASSERT(pc.first_seg_id.multi_index >= 0); + BOOST_ASSERT(pc.last_segment_index >= 0); + pc.offsetted_count = pc.last_segment_index - pc.first_seg_id.segment_index; BOOST_ASSERT(pc.offsetted_count >= 0); @@ -674,16 +898,65 @@ struct buffered_piece_collection return rob_point; } + // 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); + } + inline void calculate_robust_envelope(piece& pc) { + if (pc.offsetted_count == 0) + { + return; + } + geometry::detail::envelope::envelope_range::apply(pc.robust_ring, pc.robust_envelope); + + geometry::assign_inverse(pc.robust_offsetted_envelope); + for (int i = 0; i < pc.offsetted_count; i++) + { + geometry::expand(pc.robust_offsetted_envelope, pc.robust_ring[i]); + } + + // Take roundings into account, enlarge boxes with 1 integer + enlarge_box(pc.robust_envelope, 1); + enlarge_box(pc.robust_offsetted_envelope, 1); + } + + inline void sectionalize(piece& pc) + { + + buffered_ring<Ring> const& ring = offsetted_rings.back(); + + typedef geometry::detail::sectionalize::sectionalize_part + < + point_type, + boost::mpl::vector_c<std::size_t, 0, 1> // x,y dimension + > sectionalizer; + + // Create a ring-identifier. The source-index is the piece index + // The multi_index is as in this collection (the ring), but not used here + // The ring_index is not used + ring_identifier ring_id(pc.index, pc.first_seg_id.multi_index, -1); + + sectionalizer::apply(monotonic_sections, + boost::begin(ring) + pc.first_seg_id.segment_index, + boost::begin(ring) + pc.last_segment_index, + m_robust_policy, + ring_id, 10); } inline void finish_piece(piece& pc) { init_rescale_piece(pc, 0u); calculate_robust_envelope(pc); + sectionalize(pc); } inline void finish_piece(piece& pc, @@ -692,10 +965,16 @@ struct buffered_piece_collection const point_type& point3) { init_rescale_piece(pc, 3u); + if (pc.offsetted_count == 0) + { + return; + } + add_helper_point(pc, point1); robust_point_type mid_point = add_helper_point(pc, point2); add_helper_point(pc, point3); calculate_robust_envelope(pc); + sectionalize(pc); current_robust_ring.push_back(mid_point); } @@ -711,6 +990,7 @@ struct buffered_piece_collection robust_point_type mid_point2 = add_helper_point(pc, point2); robust_point_type mid_point1 = add_helper_point(pc, point3); add_helper_point(pc, point4); + sectionalize(pc); calculate_robust_envelope(pc); // Add mid-points in other order to current helper_ring @@ -730,10 +1010,7 @@ struct buffered_piece_collection template <typename Range> inline void add_range_to_piece(piece& pc, Range const& range, bool add_front) { - if (boost::size(range) == 0u) - { - return; - } + BOOST_ASSERT(boost::size(range) != 0u); typename Range::const_iterator it = boost::begin(range); @@ -753,10 +1030,15 @@ struct buffered_piece_collection template <typename Range> - inline void add_piece(strategy::buffer::piece_type type, Range const& range, bool decrease_segment_index_by_one) + inline void add_piece(strategy::buffer::piece_type type, Range const& range, + bool decrease_segment_index_by_one) { piece& pc = create_piece(type, decrease_segment_index_by_one); - add_range_to_piece(pc, range, offsetted_rings.back().empty()); + + if (boost::size(range) > 0u) + { + add_range_to_piece(pc, range, offsetted_rings.back().empty()); + } finish_piece(pc); } @@ -772,13 +1054,14 @@ struct buffered_piece_collection } template <typename Range> - inline void add_piece(strategy::buffer::piece_type type, point_type const& p, Range const& range) + inline void add_piece(strategy::buffer::piece_type type, + point_type const& p, Range const& range) { piece& pc = create_piece(type, true); - add_range_to_piece(pc, range, offsetted_rings.back().empty()); - if (boost::size(range) > 0) + if (boost::size(range) > 0u) { + add_range_to_piece(pc, range, offsetted_rings.back().empty()); finish_piece(pc, range.back(), p, range.front()); } else @@ -788,8 +1071,11 @@ struct buffered_piece_collection } template <typename EndcapStrategy, typename Range> - inline void add_endcap(EndcapStrategy const& strategy, Range const& range, point_type const& end_point) + inline void add_endcap(EndcapStrategy const& strategy, Range const& range, + point_type const& end_point) { + boost::ignore_unused(strategy); + if (range.empty()) { return; @@ -842,6 +1128,80 @@ struct buffered_piece_collection } } + inline bool point_coveredby_original(point_type const& point) + { + robust_point_type any_point; + geometry::recalculate(any_point, point, m_robust_policy); + + int count_in_original = 0; + + // Check of the robust point of this outputted ring is in + // any of the robust original rings + // This can go quadratic if the input has many rings, and there + // are many untouched deflated rings around + for (typename std::vector<robust_original>::const_iterator it + = robust_originals.begin(); + it != robust_originals.end(); + ++it) + { + robust_original const& original = *it; + if (detail::disjoint::disjoint_point_box(any_point, + original.m_box)) + { + continue; + } + + int const geometry_code + = detail::within::point_in_geometry(any_point, + original.m_ring); + + if (geometry_code == -1) + { + // Outside, continue + continue; + } + + // Apply for possibly nested interior rings + if (original.m_is_interior) + { + count_in_original--; + } + else if (original.m_has_interiors) + { + count_in_original++; + } + else + { + // Exterior ring without interior rings + return true; + } + } + return count_in_original > 0; + } + + // For a deflate, all rings around inner rings which are untouched + // (no intersections/turns) and which are OUTSIDE the original should + // be discarded + inline void discard_nonintersecting_deflated_rings() + { + for(typename buffered_ring_collection<buffered_ring<Ring> >::iterator it + = boost::begin(offsetted_rings); + it != boost::end(offsetted_rings); + ++it) + { + buffered_ring<Ring>& ring = *it; + if (! ring.has_intersections() + && boost::size(ring) > 0u + && geometry::area(ring) < 0) + { + if (! point_coveredby_original(geometry::range::front(ring))) + { + ring.is_untouched_outside_original = true; + } + } + } + } + inline void block_turns() { // To fix left-turn issues like #rt_u13 @@ -876,7 +1236,8 @@ struct buffered_piece_collection typedef detail::overlay::traverse < false, false, - buffered_ring_collection<buffered_ring<Ring> >, buffered_ring_collection<buffered_ring<Ring > >, + buffered_ring_collection<buffered_ring<Ring> >, + buffered_ring_collection<buffered_ring<Ring > >, backtrack_for_buffer > traverser; @@ -914,16 +1275,20 @@ struct buffered_piece_collection std::map<ring_identifier, properties> selected; - // Select all rings which do not have any self-intersection (other ones should be traversed) + // Select all rings which do not have any self-intersection + // Inner rings, for deflate, which do not have intersections, and + // which are outside originals, are skipped + // (other ones should be traversed) int index = 0; for(typename buffered_ring_collection<buffered_ring<Ring> >::const_iterator it = boost::begin(offsetted_rings); it != boost::end(offsetted_rings); ++it, ++index) { - if (! it->has_intersections()) + if (! it->has_intersections() + && ! it->is_untouched_outside_original) { ring_identifier id(0, index, -1); - selected[id] = properties(*it, true); + selected[id] = properties(*it); } } @@ -935,7 +1300,7 @@ struct buffered_piece_collection ++it, ++index) { ring_identifier id(2, index, -1); - selected[id] = properties(*it, true); + selected[id] = properties(*it); } detail::overlay::assign_parents(offsetted_rings, traversed_rings, selected, true); diff --git a/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp b/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp index 03ec598c90..9ea8bc1e85 100644 --- a/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp @@ -43,12 +43,16 @@ struct buffered_ring_collection_tag : polygonal_tag, multi_tag template <typename Ring> struct buffered_ring : public Ring { + bool has_concave; bool has_accepted_intersections; bool has_discarded_intersections; + bool is_untouched_outside_original; inline buffered_ring() - : has_accepted_intersections(false) + : has_concave(false) + , has_accepted_intersections(false) , has_discarded_intersections(false) + , is_untouched_outside_original(false) {} inline bool discarded() const diff --git a/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp b/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp index 395921ccaa..6a3daa282e 100644 --- a/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp +++ b/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp @@ -16,6 +16,7 @@ #include <boost/geometry/algorithms/detail/disjoint/box_box.hpp> #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp> #include <boost/geometry/algorithms/detail/overlay/get_turn_info.hpp> +#include <boost/geometry/algorithms/detail/sections/section_functions.hpp> namespace boost { namespace geometry @@ -27,32 +28,16 @@ namespace detail { namespace buffer { -struct piece_get_box -{ - template <typename Box, typename Piece> - static inline void apply(Box& total, Piece const& piece) - { - geometry::expand(total, piece.robust_envelope); - } -}; - -struct piece_ovelaps_box -{ - template <typename Box, typename Piece> - static inline bool apply(Box const& box, Piece const& piece) - { - return ! geometry::detail::disjoint::disjoint_box_box(box, piece.robust_envelope); - } -}; - template < + typename Pieces, typename Rings, typename Turns, typename RobustPolicy > class piece_turn_visitor { + Pieces const& m_pieces; Rings const& m_rings; Turns& m_turns; RobustPolicy const& m_robust_policy; @@ -69,6 +54,17 @@ class piece_turn_visitor || piece1.index == piece2.right_index; } + template <typename Piece> + inline bool is_on_same_convex_ring(Piece const& piece1, Piece const& piece2) const + { + if (piece1.first_seg_id.multi_index != piece2.first_seg_id.multi_index) + { + return false; + } + + return ! m_rings[piece1.first_seg_id.multi_index].has_concave; + } + template <typename Range, typename Iterator> inline void move_to_next_point(Range const& range, Iterator& next) const { @@ -93,46 +89,110 @@ class piece_turn_visitor return result; } - template <typename Piece> - inline void calculate_turns(Piece const& piece1, Piece const& piece2) + template <std::size_t Dimension, typename Iterator, typename Box> + inline void move_begin_iterator(Iterator& it_begin, Iterator it_beyond, + int& index, int dir, Box const& other_bounding_box) + { + for(; it_begin != it_beyond + && it_begin + 1 != it_beyond + && detail::section::preceding<Dimension>(dir, *(it_begin + 1), + other_bounding_box, m_robust_policy); + ++it_begin, index++) + {} + } + + template <std::size_t Dimension, typename Iterator, typename Box> + inline void move_end_iterator(Iterator it_begin, Iterator& it_beyond, + int dir, Box const& other_bounding_box) + { + while (it_beyond != it_begin + && it_beyond - 1 != it_begin + && it_beyond - 2 != it_begin) + { + if (detail::section::exceeding<Dimension>(dir, *(it_beyond - 2), + other_bounding_box, m_robust_policy)) + { + --it_beyond; + } + else + { + return; + } + } + } + + template <typename Piece, typename Section> + inline void calculate_turns(Piece const& piece1, Piece const& piece2, + Section const& section1, Section const& section2) { typedef typename boost::range_value<Rings const>::type ring_type; typedef typename boost::range_value<Turns const>::type turn_type; typedef typename boost::range_iterator<ring_type const>::type iterator; - segment_identifier seg_id1 = piece1.first_seg_id; - segment_identifier seg_id2 = piece2.first_seg_id; - - if (seg_id1.segment_index < 0 || seg_id2.segment_index < 0) + int const piece1_first_index = piece1.first_seg_id.segment_index; + int const piece2_first_index = piece2.first_seg_id.segment_index; + if (piece1_first_index < 0 || piece2_first_index < 0) { return; } - ring_type const& ring1 = m_rings[seg_id1.multi_index]; - iterator it1_first = boost::begin(ring1) + seg_id1.segment_index; - iterator it1_last = boost::begin(ring1) + piece1.last_segment_index; - - ring_type const& ring2 = m_rings[seg_id2.multi_index]; - iterator it2_first = boost::begin(ring2) + seg_id2.segment_index; - iterator it2_last = boost::begin(ring2) + piece2.last_segment_index; + // Get indices of part of offsetted_rings for this monotonic section: + int const sec1_first_index = piece1_first_index + section1.begin_index; + int const sec2_first_index = piece2_first_index + section2.begin_index; + + // index of last point in section, beyond-end is one further + int const sec1_last_index = piece1_first_index + section1.end_index; + int const sec2_last_index = piece2_first_index + section2.end_index; + + // get geometry and iterators over these sections + ring_type const& ring1 = m_rings[piece1.first_seg_id.multi_index]; + iterator it1_first = boost::begin(ring1) + sec1_first_index; + iterator it1_beyond = boost::begin(ring1) + sec1_last_index + 1; + + ring_type const& ring2 = m_rings[piece2.first_seg_id.multi_index]; + iterator it2_first = boost::begin(ring2) + sec2_first_index; + iterator it2_beyond = boost::begin(ring2) + sec2_last_index + 1; + + // Set begin/end of monotonic ranges, in both x/y directions + int index1 = sec1_first_index; + move_begin_iterator<0>(it1_first, it1_beyond, index1, + section1.directions[0], section2.bounding_box); + move_end_iterator<0>(it1_first, it1_beyond, + section1.directions[0], section2.bounding_box); + move_begin_iterator<1>(it1_first, it1_beyond, index1, + section1.directions[1], section2.bounding_box); + move_end_iterator<1>(it1_first, it1_beyond, + section1.directions[1], section2.bounding_box); + + int index2 = sec2_first_index; + move_begin_iterator<0>(it2_first, it2_beyond, index2, + section2.directions[0], section1.bounding_box); + move_end_iterator<0>(it2_first, it2_beyond, + section2.directions[0], section1.bounding_box); + move_begin_iterator<1>(it2_first, it2_beyond, index2, + section2.directions[1], section1.bounding_box); + move_end_iterator<1>(it2_first, it2_beyond, + section2.directions[1], section1.bounding_box); turn_type the_model; the_model.operations[0].piece_index = piece1.index; the_model.operations[0].seg_id = piece1.first_seg_id; + the_model.operations[0].seg_id.segment_index = index1; // override iterator it1 = it1_first; for (iterator prev1 = it1++; - it1 != it1_last; + it1 != it1_beyond; prev1 = it1++, the_model.operations[0].seg_id.segment_index++) { the_model.operations[1].piece_index = piece2.index; the_model.operations[1].seg_id = piece2.first_seg_id; + the_model.operations[1].seg_id.segment_index = index2; // override iterator next1 = next_point(ring1, it1); iterator it2 = it2_first; for (iterator prev2 = it2++; - it2 != it2_last; + it2 != it2_beyond; prev2 = it2++, the_model.operations[1].seg_id.segment_index++) { iterator next2 = next_point(ring2, it2); @@ -158,26 +218,36 @@ class piece_turn_visitor public: - piece_turn_visitor(Rings const& ring_collection, + piece_turn_visitor(Pieces const& pieces, + Rings const& ring_collection, Turns& turns, RobustPolicy const& robust_policy) - : m_rings(ring_collection) + : m_pieces(pieces) + , m_rings(ring_collection) , m_turns(turns) , m_robust_policy(robust_policy) {} - template <typename Piece> - inline void apply(Piece const& piece1, Piece const& piece2, + template <typename Section> + inline void apply(Section const& section1, Section const& section2, bool first = true) { boost::ignore_unused_variable_warning(first); - if ( is_adjacent(piece1, piece2) - || detail::disjoint::disjoint_box_box(piece1.robust_envelope, - piece2.robust_envelope)) + + typedef typename boost::range_value<Pieces const>::type piece_type; + piece_type const& piece1 = m_pieces[section1.ring_id.source_index]; + piece_type const& piece2 = m_pieces[section2.ring_id.source_index]; + + if ( piece1.index == piece2.index + || is_adjacent(piece1, piece2) + || is_on_same_convex_ring(piece1, piece2) + || detail::disjoint::disjoint_box_box(section1.bounding_box, + section2.bounding_box) ) { return; } - calculate_turns(piece1, piece2); + + calculate_turns(piece1, piece2, section1, section2); } }; diff --git a/boost/geometry/algorithms/detail/buffer/turn_in_input.hpp b/boost/geometry/algorithms/detail/buffer/turn_in_input.hpp deleted file mode 100644 index 2b1c33d291..0000000000 --- a/boost/geometry/algorithms/detail/buffer/turn_in_input.hpp +++ /dev/null @@ -1,98 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) - -// Copyright (c) 2012-2014 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_BUFFER_TURN_IN_INPUT_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_INPUT_HPP - -#include <boost/geometry/core/tags.hpp> -#include <boost/geometry/algorithms/covered_by.hpp> - - -namespace boost { namespace geometry -{ - -#ifndef DOXYGEN_NO_DETAIL -namespace detail { namespace buffer -{ - -// Checks if an turn/intersection point is inside (or covered by) the input geometry - -template <typename Tag> -struct turn_in_input -{ -}; - -template <> -struct turn_in_input<polygon_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& point, Geometry const& geometry) - { - return geometry::covered_by(point, geometry) ? 1 : -1; - } -}; - -template <> -struct turn_in_input<linestring_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& , Geometry const& ) - { - return 0; - } -}; - -template <> -struct turn_in_input<point_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& , Geometry const& ) - { - return 0; - } -}; - -template <> -struct turn_in_input<multi_polygon_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& point, Geometry const& geometry) - { - return geometry::covered_by(point, geometry) ? 1 : -1; - } -}; - -template <> -struct turn_in_input<multi_linestring_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& , Geometry const& ) - { - return 0; - } -}; - -template <> -struct turn_in_input<multi_point_tag> -{ - template <typename Point, typename Geometry> - static inline int apply(Point const& , Geometry const& ) - { - return 0; - } -}; - - -}} // namespace detail::buffer -#endif // DOXYGEN_NO_DETAIL - - - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_INPUT_HPP diff --git a/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp b/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp new file mode 100644 index 0000000000..05fc6df8ff --- /dev/null +++ b/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp @@ -0,0 +1,268 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2014 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_BUFFER_TURN_IN_ORIGINAL_VISITOR +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_ORIGINAL_VISITOR + + +#include <boost/core/ignore_unused.hpp> + +#include <boost/geometry/algorithms/expand.hpp> +#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> +#include <boost/geometry/strategies/buffer.hpp> + + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace buffer +{ + +struct original_get_box +{ + template <typename Box, typename Original> + static inline void apply(Box& total, Original const& original) + { + geometry::expand(total, original.m_box); + } +}; + +struct original_ovelaps_box +{ + template <typename Box, typename Original> + static inline bool apply(Box const& box, Original const& original) + { + return ! detail::disjoint::disjoint_box_box(box, original.m_box); + } +}; + +struct include_turn_policy +{ + template <typename Turn> + static inline bool apply(Turn const& turn) + { + return turn.location == location_ok; + } +}; + +struct turn_in_original_ovelaps_box +{ + template <typename Box, typename Turn> + static inline bool apply(Box const& box, Turn const& turn) + { + if (turn.location != location_ok || turn.within_original) + { + // Skip all points already processed + return false; + } + + return ! geometry::detail::disjoint::disjoint_point_box( + turn.robust_point, box); + } +}; + +//! Check if specified is in range of specified iterators +//! Return value of strategy (true if we can bail out) +template +< + typename Strategy, + typename State, + typename Point, + typename Iterator +> +inline bool point_in_range(Strategy& strategy, State& state, + Point const& point, Iterator begin, Iterator end) +{ + boost::ignore_unused(strategy); + + Iterator it = begin; + for (Iterator previous = it++; it != end; ++previous, ++it) + { + if (! strategy.apply(point, *previous, *it, state)) + { + // We're probably on the boundary + return false; + } + } + return true; +} + +template +< + typename Strategy, + typename State, + typename Point, + typename CoordinateType, + typename Iterator +> +inline bool point_in_section(Strategy& strategy, State& state, + Point const& point, CoordinateType const& point_y, + Iterator begin, Iterator end, + int direction) +{ + if (direction == 0) + { + // Not a monotonic section, or no change in Y-direction + return point_in_range(strategy, state, point, begin, end); + } + + // We're in a monotonic section in y-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); + + if (direction == 1 && point_y < previous_y) + { + // Section goes upwards, y increases, point is is below section + return true; + } + else if (direction == -1 && point_y > previous_y) + { + // Section goes downwards, y decreases, point is above section + return true; + } + + if (! strategy.apply(point, *previous, *it, state)) + { + // We're probably on the boundary + return false; + } + } + return true; +} + + +template <typename Point, typename Original> +inline int point_in_original(Point const& point, Original const& original) +{ + typedef strategy::within::winding<Point> strategy_type; + + typename strategy_type::state_type state; + strategy_type strategy; + + if (boost::size(original.m_sections) == 0 + || boost::size(original.m_ring) - boost::size(original.m_sections) < 16) + { + // There are no sections, or it does not profit to walk over sections + // instead of over points. Boundary of 16 is arbitrary but can influence + // performance + point_in_range(strategy, state, point, + original.m_ring.begin(), original.m_ring.end()); + return strategy.result(state); + } + + typedef typename Original::sections_type sections_type; + typedef typename boost::range_iterator<sections_type const>::type iterator_type; + 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); + + // Walk through all monotonic sections of this original + for (iterator_type it = boost::begin(original.m_sections); + it != boost::end(original.m_sections); + ++it) + { + section_type const& section = *it; + + 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)) + { + // y-coordinate of point overlaps with section + if (! point_in_section(strategy, state, point, point_y, + boost::begin(original.m_ring) + section.begin_index, + boost::begin(original.m_ring) + section.end_index + 1, + section.directions[0])) + { + // We're probably on the boundary + break; + } + } + } + + return strategy.result(state); +} + + +template <typename Turns> +class turn_in_original_visitor +{ +public: + turn_in_original_visitor(Turns& turns) + : m_mutable_turns(turns) + {} + + template <typename Turn, typename Original> + inline void apply(Turn const& turn, Original const& original, bool first = true) + { + boost::ignore_unused_variable_warning(first); + + if (turn.location != location_ok || turn.within_original) + { + // Skip all points already processed + return; + } + + if (geometry::disjoint(turn.robust_point, original.m_box)) + { + // Skip all disjoint + return; + } + + int const code = point_in_original(turn.robust_point, original); + + if (code == -1) + { + return; + } + + Turn& mutable_turn = m_mutable_turns[turn.turn_index]; + + if (code == 0) + { + // On border of original: always discard + mutable_turn.location = location_discard; + } + + // Point is inside an original ring + if (original.m_is_interior) + { + mutable_turn.count_in_original--; + } + else if (original.m_has_interiors) + { + mutable_turn.count_in_original++; + } + else + { + // It is an exterior ring and there are no interior rings. + // Then we are completely ready with this turn + mutable_turn.within_original = true; + mutable_turn.count_in_original = 1; + } + } + +private : + Turns& m_mutable_turns; +}; + + +}} // namespace detail::buffer +#endif // DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_ORIGINAL_VISITOR 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 c02f56de3f..8803efdec9 100644 --- a/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp +++ b/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp @@ -9,19 +9,23 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_PIECE_VISITOR #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_TURN_IN_PIECE_VISITOR + +#include <boost/core/ignore_unused.hpp> + #include <boost/range.hpp> #include <boost/geometry/arithmetic/dot_product.hpp> +#include <boost/geometry/algorithms/assign.hpp> +#include <boost/geometry/algorithms/comparable_distance.hpp> #include <boost/geometry/algorithms/equals.hpp> #include <boost/geometry/algorithms/expand.hpp> #include <boost/geometry/algorithms/detail/disjoint/point_box.hpp> +#include <boost/geometry/algorithms/detail/disjoint/box_box.hpp> #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp> #include <boost/geometry/algorithms/detail/overlay/get_turn_info.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/strategies/buffer.hpp> -#include <boost/geometry/geometries/linestring.hpp> -#include <boost/geometry/algorithms/comparable_distance.hpp> namespace boost { namespace geometry { @@ -31,6 +35,33 @@ namespace boost { namespace geometry namespace detail { namespace buffer { +struct piece_get_box +{ + template <typename Box, typename Piece> + static inline void apply(Box& total, Piece const& piece) + { + geometry::expand(total, piece.robust_envelope); + } +}; + +struct piece_ovelaps_box +{ + template <typename Box, typename Piece> + static inline bool apply(Box const& box, Piece const& piece) + { + if (piece.type == strategy::buffer::buffered_flat_end + || piece.type == strategy::buffer::buffered_concave) + { + // Turns cannot be inside a flat end (though they can be on border) + // Neither we need to check if they are inside concave helper pieces + + // Skip all pieces not used as soon as possible + return false; + } + + return ! geometry::detail::disjoint::disjoint_box_box(box, piece.robust_envelope); + } +}; struct turn_get_box { @@ -46,96 +77,422 @@ struct turn_ovelaps_box template <typename Box, typename Turn> static inline bool apply(Box const& box, Turn const& turn) { - return ! dispatch::disjoint - < - typename Turn::robust_point_type, - Box - >::apply(turn.robust_point, box); + return ! geometry::detail::disjoint::disjoint_point_box(turn.robust_point, box); } }; -template <typename Turns, typename Pieces> -class turn_in_piece_visitor + +enum analyse_result { - Turns& m_turns; // because partition is currently operating on const input only - Pieces const& m_pieces; // to check for piece-type + analyse_unknown, + analyse_continue, + analyse_disjoint, + analyse_within, + analyse_on_original_boundary, + analyse_on_offsetted, + analyse_near_offsetted +}; - typedef boost::long_long_type calculation_type; +template <typename Point> +inline bool in_box(Point const& previous, + Point const& current, Point const& point) +{ + // Get its box (TODO: this can be prepared-on-demand later) + typedef geometry::model::box<Point> box_type; + box_type box; + geometry::assign_inverse(box); + geometry::expand(box, previous); + geometry::expand(box, current); + + return geometry::covered_by(point, box); +} + +template <typename Point, typename Turn> +inline analyse_result check_segment(Point const& previous, + Point const& current, Turn const& turn, + bool from_monotonic) +{ + typedef typename strategy::side::services::default_strategy + < + typename cs_tag<Point>::type + >::type side_strategy; + typedef typename geometry::coordinate_type<Point>::type coordinate_type; + + coordinate_type const twice_area + = side_strategy::template side_value + < + coordinate_type, + coordinate_type + >(previous, current, turn.robust_point); - template <typename Point> - static inline bool projection_on_segment(Point const& subject, Point const& p, Point const& q) + if (twice_area == 0) { - typedef Point vector_type; - typedef typename geometry::coordinate_type<Point>::type coordinate_type; + // Collinear, only on segment if it is covered by its bbox + if (in_box(previous, current, turn.robust_point)) + { + return analyse_on_offsetted; + } + } + else if (twice_area < 0) + { + // It is in the triangle right-of the segment where the + // segment is the hypothenusa. Check if it is close + // (within rounding-area) + if (twice_area * twice_area < geometry::comparable_distance(previous, current) + && in_box(previous, current, turn.robust_point)) + { + return analyse_near_offsetted; + } + else if (from_monotonic) + { + return analyse_within; + } + } + else if (twice_area > 0 && from_monotonic) + { + // Left of segment + return analyse_disjoint; + } + + // Not monotonic, on left or right side: continue analysing + return analyse_continue; +} - vector_type v = q; - vector_type w = subject; - subtract_point(v, p); - subtract_point(w, p); - coordinate_type const zero = coordinate_type(); - coordinate_type const c1 = dot_product(w, v); +class analyse_turn_wrt_point_piece +{ +public : + template <typename Turn, typename Piece> + static inline analyse_result apply(Turn const& turn, Piece const& piece) + { + typedef typename Piece::section_type section_type; + typedef typename Turn::robust_point_type point_type; + typedef typename geometry::coordinate_type<point_type>::type coordinate_type; - if (c1 < zero) + coordinate_type const point_y = geometry::get<1>(turn.robust_point); + + typedef strategy::within::winding<point_type> strategy_type; + + typename strategy_type::state_type state; + strategy_type strategy; + boost::ignore_unused(strategy); + + for (std::size_t s = 0; s < piece.sections.size(); s++) { - return false; + section_type const& section = piece.sections[s]; + // If point within vertical 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) + { + for (int i = section.begin_index + 1; i <= section.end_index; i++) + { + point_type const& previous = piece.robust_ring[i - 1]; + point_type const& current = piece.robust_ring[i]; + + analyse_result code = check_segment(previous, current, turn, false); + if (code != analyse_continue) + { + return code; + } + + // Get the state (to determine it is within), we don't have + // to cover the on-segment case (covered above) + strategy.apply(turn.robust_point, previous, current, state); + } + } } - coordinate_type const c2 = dot_product(v, v); - if (c2 < c1) + + int const code = strategy.result(state); + if (code == 1) { - return false; + return analyse_within; + } + else if (code == -1) + { + return analyse_disjoint; } - return true; + // Should normally not occur - on-segment is covered + return analyse_unknown; } - template <typename Point, typename Piece> - inline bool on_offsetted(Point const& point, Piece const& piece) const +}; + +class analyse_turn_wrt_piece +{ + template <typename Point, typename Turn> + static inline analyse_result check_helper_segment(Point const& s1, + Point const& s2, Turn const& turn, + bool is_original, + Point const& offsetted) { typedef typename strategy::side::services::default_strategy < typename cs_tag<Point>::type >::type side_strategy; - geometry::equal_to<Point> comparator; + + switch(side_strategy::apply(s1, s2, turn.robust_point)) + { + case 1 : + return analyse_disjoint; // left of segment + case 0 : + { + // If is collinear, either on segment or before/after + typedef geometry::model::box<Point> box_type; + + box_type box; + geometry::assign_inverse(box); + geometry::expand(box, s1); + geometry::expand(box, s2); + + if (geometry::covered_by(turn.robust_point, box)) + { + // It is on the segment + if (! is_original + && geometry::comparable_distance(turn.robust_point, offsetted) <= 1) + { + // It is close to the offsetted-boundary, take + // any rounding-issues into account + return analyse_near_offsetted; + } + + // Points on helper-segments are considered as within + // Points on original boundary are processed differently + return is_original + ? analyse_on_original_boundary + : analyse_within; + } + + // It is collinear but not on the segment. Because these + // segments are convex, it is outside + // Unless the offsetted ring is collinear or concave w.r.t. + // helper-segment but that scenario is not yet supported + return analyse_disjoint; + } + break; + } + + // right of segment + return analyse_continue; + } + + template <typename Turn, typename Piece> + static inline analyse_result check_helper_segments(Turn const& turn, Piece const& piece) + { + typedef typename Turn::robust_point_type point_type; + geometry::equal_to<point_type> comparator; + + point_type points[4]; + + int helper_count = piece.robust_ring.size() - piece.offsetted_count; + if (helper_count == 4) + { + for (int i = 0; i < 4; i++) + { + points[i] = piece.robust_ring[piece.offsetted_count + i]; + } + } + else if (helper_count == 3) + { + // Triangular piece, assign points but assign second twice + for (int i = 0; i < 4; i++) + { + int index = i < 2 ? i : i - 1; + points[i] = piece.robust_ring[piece.offsetted_count + index]; + } + } + else + { + // Some pieces (e.g. around points) do not have helper segments. + // Others should have 3 (join) or 4 (side) + return analyse_continue; + } + + // First check point-equality + point_type const& point = turn.robust_point; + if (comparator(point, points[0]) || comparator(point, points[3])) + { + return analyse_on_offsetted; + } + if (comparator(point, points[1]) || comparator(point, points[2])) + { + return analyse_on_original_boundary; + } + + // Right side of the piece + analyse_result result + = check_helper_segment(points[0], points[1], turn, + false, points[0]); + if (result != analyse_continue) + { + return result; + } + + // Left side of the piece + result = check_helper_segment(points[2], points[3], turn, + false, points[3]); + if (result != analyse_continue) + { + return result; + } + + if (! comparator(points[1], points[2])) + { + // Side of the piece at side of original geometry + result = check_helper_segment(points[1], points[2], turn, + true, point); + if (result != analyse_continue) + { + return result; + } + } + + // We are within the \/ or |_| shaped piece, where the top is the + // offsetted ring. + if (! geometry::covered_by(point, piece.robust_offsetted_envelope)) + { + // Not in offsetted-area. This makes a cheap check possible + typedef typename strategy::side::services::default_strategy + < + typename cs_tag<point_type>::type + >::type side_strategy; + + switch(side_strategy::apply(points[3], points[0], point)) + { + case 1 : return analyse_disjoint; + case -1 : return analyse_within; + case 0 : return analyse_disjoint; + } + } + + return analyse_continue; + } + + template <typename Turn, typename Piece, typename Compare> + static inline analyse_result check_monotonic(Turn const& turn, Piece const& piece, Compare const& compare) + { + typedef typename Piece::piece_robust_ring_type ring_type; + typedef typename ring_type::const_iterator it_type; + it_type end = piece.robust_ring.begin() + piece.offsetted_count; + it_type it = std::lower_bound(piece.robust_ring.begin(), + end, + turn.robust_point, + compare); + + if (it != end + && it != piece.robust_ring.begin()) + { + // iterator points to point larger than point + // w.r.t. specified direction, and prev points to a point smaller + // We now know if it is inside/outside + it_type prev = it - 1; + return check_segment(*prev, *it, turn, true); + } + return analyse_continue; + } + +public : + template <typename Turn, typename Piece> + static inline analyse_result apply(Turn const& turn, Piece const& piece) + { + typedef typename Turn::robust_point_type point_type; + analyse_result code = check_helper_segments(turn, piece); + if (code != analyse_continue) + { + return code; + } + + geometry::equal_to<point_type> comparator; + + if (piece.offsetted_count > 8) + { + // If the offset contains some points and is monotonic, we try + // to avoid walking all points linearly. + // We try it only once. + if (piece.is_monotonic_increasing[0]) + { + code = check_monotonic(turn, piece, geometry::less<point_type, 0>()); + if (code != analyse_continue) return code; + } + else if (piece.is_monotonic_increasing[1]) + { + code = check_monotonic(turn, piece, geometry::less<point_type, 1>()); + if (code != analyse_continue) return code; + } + else if (piece.is_monotonic_decreasing[0]) + { + code = check_monotonic(turn, piece, geometry::greater<point_type, 0>()); + if (code != analyse_continue) return code; + } + else if (piece.is_monotonic_decreasing[1]) + { + code = check_monotonic(turn, piece, geometry::greater<point_type, 1>()); + if (code != analyse_continue) return code; + } + } + + // It is small or not monotonic, walk linearly through offset + // TODO: this will be combined with winding strategy for (int i = 1; i < piece.offsetted_count; i++) { - Point const& previous = piece.robust_ring[i - 1]; - Point const& current = piece.robust_ring[i]; + point_type const& previous = piece.robust_ring[i - 1]; + point_type const& current = piece.robust_ring[i]; - // The robust ring contains duplicates, avoid applying side on them (will be 0) + // The robust ring can contain duplicates + // (on which any side or side-value would return 0) if (! comparator(previous, current)) { - int const side = side_strategy::apply(previous, current, point); - if (side == 0) + code = check_segment(previous, current, turn, false); + if (code != analyse_continue) { - // Collinear, check if projection falls on it - if (projection_on_segment(point, previous, current)) - { - return true; - } + return code; } } } - return false; + + return analyse_unknown; } - template <typename Point, typename Piece> - static inline - calculation_type comparable_distance_from_offsetted(Point const& point, - Piece const& piece) +}; + + +template <typename Turns, typename Pieces> +class turn_in_piece_visitor +{ + Turns& m_turns; // because partition is currently operating on const input only + Pieces const& m_pieces; // to check for piece-type + + template <typename Operation, typename Piece> + inline bool skip(Operation const& op, Piece const& piece) const { - // TODO: pass subrange to dispatch to avoid making copy - geometry::model::linestring<Point> ls; - std::copy(piece.robust_ring.begin(), - piece.robust_ring.begin() + piece.offsetted_count, - std::back_inserter(ls)); - typename default_comparable_distance_result<Point, Point>::type - const comp = geometry::comparable_distance(point, ls); - - return static_cast<calculation_type>(comp); + if (op.piece_index == piece.index) + { + return true; + } + Piece const& pc = m_pieces[op.piece_index]; + if (pc.left_index == piece.index || pc.right_index == piece.index) + { + if (pc.type == strategy::buffer::buffered_flat_end) + { + // If it is a flat end, don't compare against its neighbor: + // it will always be located on one of the helper segments + return true; + } + if (pc.type == strategy::buffer::buffered_concave) + { + // If it is concave, the same applies: the IP will be + // located on one of the helper segments + return true; + } + } + + return false; } + public: inline turn_in_piece_visitor(Turns& turns, Pieces const& pieces) @@ -157,81 +514,55 @@ public: if (piece.type == strategy::buffer::buffered_flat_end || piece.type == strategy::buffer::buffered_concave) { - // Turns cannot be inside a flat end (though they can be on border) - // Neither we need to check if they are inside concave helper pieces + // Turns cannot be located within flat-end or concave pieces return; } - bool neighbour = false; - for (int i = 0; i < 2; i++) - { - // Don't compare against one of the two source-pieces - if (turn.operations[i].piece_index == piece.index) - { - return; - } - - typename boost::range_value<Pieces>::type const& pc - = m_pieces[turn.operations[i].piece_index]; - if (pc.left_index == piece.index - || pc.right_index == piece.index) - { - if (pc.type == strategy::buffer::buffered_flat_end) - { - // If it is a flat end, don't compare against its neighbor: - // it will always be located on one of the helper segments - return; - } - neighbour = true; - } - } - - int geometry_code = detail::within::point_in_geometry(turn.robust_point, piece.robust_ring); - - if (geometry_code == -1) + if (! geometry::covered_by(turn.robust_point, piece.robust_envelope)) { + // Easy check: if the turn is not in the envelope, we can safely return return; } - if (geometry_code == 0 && neighbour) + + if (skip(turn.operations[0], piece) || skip(turn.operations[1], piece)) { return; } + // TODO: mutable_piece to make some on-demand preparations in analyse + analyse_result analyse_code = + piece.type == geometry::strategy::buffer::buffered_point + ? analyse_turn_wrt_point_piece::apply(turn, piece) + : analyse_turn_wrt_piece::apply(turn, piece); + Turn& mutable_turn = m_turns[turn.turn_index]; - if (geometry_code == 0) + switch(analyse_code) { - // If it is on the border and they are not neighbours, it should be - // on the offsetted ring - - if (! on_offsetted(turn.robust_point, piece)) - { - // It is on the border but not on the offsetted ring. - // Then it is somewhere on the helper-segments - // Classify it as "within" - geometry_code = 1; - mutable_turn.count_on_helper++; // can still become "near_offsetted" - } - else - { + case analyse_disjoint : + return; + case analyse_on_offsetted : mutable_turn.count_on_offsetted++; // value is not used anymore - } + return; + case analyse_on_original_boundary : + mutable_turn.count_on_original_boundary++; + return; + case analyse_within : + mutable_turn.count_within++; + return; + case analyse_near_offsetted : + mutable_turn.count_within_near_offsetted++; + return; + default : + break; } + // TODO: this point_in_geometry is a performance-bottleneck here and + // will be replaced completely by extending analyse_piece functionality + int geometry_code = detail::within::point_in_geometry(turn.robust_point, piece.robust_ring); + if (geometry_code == 1) { - calculation_type const distance - = comparable_distance_from_offsetted(turn.robust_point, piece); - if (distance >= 4) - { - // This is too far from the border, it counts as "really within" - mutable_turn.count_within++; - } - else - { - // Other points count as still "on border" because they might be - // travelled through, but not used as starting point - mutable_turn.count_within_near_offsetted++; - } + mutable_turn.count_within++; } } }; |