From 447fd7edd2d94872bf3de63610ccb25d8914b205 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Mon, 9 Jun 2014 15:07:47 +0200 Subject: [PATCH] [buffer] get occupation vectors by rescaled offsetted points This fixes the last case rt_p20 Check if the occupation vectors are short (length 1). This indicates a rounding issue. If so, map again but use neighbouring cells. Alas we have to do this, but still better than the former FP implementation. Also, we map only points on offsetted borders now, and return if the map is empty, to improve performance. --- .../buffer/multi_polygon_buffer.cpp | 2 - .../algorithms/detail/get_left_turns.hpp | 39 +++-- .../algorithms/detail/occupation_info.hpp | 123 ++++++-------- .../algorithms/buffer/buffer_policies.hpp | 5 + .../buffer/buffered_piece_collection.hpp | 159 ++++++++++++++---- .../buffered_piece_collection_with_mapper.hpp | 7 +- 6 files changed, 210 insertions(+), 125 deletions(-) diff --git a/extensions/test/algorithms/buffer/multi_polygon_buffer.cpp b/extensions/test/algorithms/buffer/multi_polygon_buffer.cpp index 554dfdc7e..f75347396 100644 --- a/extensions/test/algorithms/buffer/multi_polygon_buffer.cpp +++ b/extensions/test/algorithms/buffer/multi_polygon_buffer.cpp @@ -335,9 +335,7 @@ void test_all() test_one("rt_p17", rt_p17, 25.3137, 1.0); test_one("rt_p18", rt_p18, 23.3137, 1.0); test_one("rt_p19", rt_p19, 25.5637, 1.0); -#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS) test_one("rt_p20", rt_p20, 25.4853, 1.0); -#endif test_one("rt_p21", rt_p21, 17.1716, 1.0); test_one("rt_p22", rt_p22, 26.5711, 1.0); diff --git a/include/boost/geometry/algorithms/detail/get_left_turns.hpp b/include/boost/geometry/algorithms/detail/get_left_turns.hpp index 527039fbe..45f5f66ea 100644 --- a/include/boost/geometry/algorithms/detail/get_left_turns.hpp +++ b/include/boost/geometry/algorithms/detail/get_left_turns.hpp @@ -79,6 +79,14 @@ inline int get_quadrant(Vector const& vector) ; } +template +inline int squared_length(Vector const& vector) +{ + return geometry::get<0>(vector) * geometry::get<0>(vector) + + geometry::get<1>(vector) * geometry::get<1>(vector) + ; +} + template struct angle_less @@ -93,13 +101,6 @@ struct angle_less : m_origin(origin) {} - inline int length(vector_type const& v) const - { - return geometry::get<0>(v) * geometry::get<0>(v) - + geometry::get<1>(v) * geometry::get<1>(v) - ; - } - template inline bool operator()(Angle const& p, Angle const& q) const { @@ -128,11 +129,11 @@ struct angle_less return int(p.incoming) < int(q.incoming); } // Same quadrant/side/direction, return longest first - int const length_p = length(pv); - int const length_q = length(qv); + int const length_p = squared_length(pv); + int const length_q = squared_length(qv); if (length_p != length_q) { - return length(pv) > length(qv); + return squared_length(pv) > squared_length(qv); } // They are still the same. Just compare on seg_id return p.seg_id < q.seg_id; @@ -259,6 +260,24 @@ inline void block_turns_on_right_sides(AngleTurnCollection const& turns, } } +template +inline bool has_rounding_issues(AngleCollection const& angles, Point const& origin) +{ + for (typename boost::range_iterator::type it = + angles.begin(); it != angles.end(); ++it) + { + // Vector origin -> p and origin -> q + typedef Point vector_type; + vector_type v = it->point; + geometry::subtract_point(v, origin); + return geometry::math::abs(geometry::get<0>(v)) <= 1 + || geometry::math::abs(geometry::get<1>(v)) <= 1 + ; + } + return false; +} + + } // namespace left_turns } // namespace detail diff --git a/include/boost/geometry/algorithms/detail/occupation_info.hpp b/include/boost/geometry/algorithms/detail/occupation_info.hpp index 440a34fb0..7980c6d7b 100644 --- a/include/boost/geometry/algorithms/detail/occupation_info.hpp +++ b/include/boost/geometry/algorithms/detail/occupation_info.hpp @@ -29,36 +29,6 @@ namespace boost { namespace geometry namespace detail { -template -< - typename Iterator, - typename Vector, - typename RobustPolicy, - typename RobustPoint -> -inline Iterator advance_circular(Iterator it, Vector const& vector, - RobustPolicy const& robust_policy, - RobustPoint& robust_point, - segment_identifier& seg_id, - bool forward = true) -{ - int const increment = forward ? 1 : -1; - if (it == boost::begin(vector) && increment < 0) - { - it = boost::end(vector); - seg_id.segment_index = boost::size(vector); - } - it += increment; - seg_id.segment_index += increment; - if (it == boost::end(vector)) - { - seg_id.segment_index = 0; - it = boost::begin(vector); - } - geometry::recalculate(robust_point, *it, robust_policy); - return it; -} - template struct angle_info { @@ -146,70 +116,77 @@ public : detail::left_turns::get_left_turns(angles, origin, turns_to_keep); } + + template + inline bool has_rounding_issues(RobustPoint const& origin) const + { + return detail::left_turns::has_rounding_issues(angles, origin); + } }; +template +inline void move_index(Pieces const& pieces, int& index, int& piece_index, int direction) +{ + BOOST_ASSERT(direction == 1 || direction == -1); + BOOST_ASSERT(piece_index >= 0 && piece_index < boost::size(pieces)); + BOOST_ASSERT(index >= 0 && index < boost::size(pieces[piece_index].robust_ring)); + + index += direction; + if (direction == -1 && index < 0) + { + piece_index--; + if (piece_index < 0) + { + piece_index = boost::size(pieces) - 1; + } + index = boost::size(pieces[piece_index].robust_ring) - 1; + } + if (direction == 1 && index >= boost::size(pieces[piece_index].robust_ring)) + { + piece_index++; + if (piece_index >= boost::size(pieces)) + { + piece_index = 0; + } + index = 0; + } +} + template < typename RobustPoint, - typename RobustPolicy, - typename Ring, + typename Turn, + typename Pieces, typename Info > inline void add_incoming_and_outgoing_angles( RobustPoint const& intersection_point, // rescaled - Ring const& ring, // non-rescaled - RobustPolicy const& robust_policy, + Turn const& turn, + Pieces const& pieces, // using rescaled offsets of it int turn_index, int operation_index, segment_identifier seg_id, Info& info) { - typedef typename boost::range_iterator - < - Ring const - >::type iterator_type; - - int const n = boost::size(ring); - if (seg_id.segment_index >= n || seg_id.segment_index < 0) - { - return; - } - segment_identifier real_seg_id = seg_id; - iterator_type it = boost::begin(ring) + seg_id.segment_index; - RobustPoint current; - geometry::recalculate(current, *it, robust_policy); - - // TODO: we can add this points in get_turn_info/assign already - // TODO: if we use turn-info ("to", "middle"), we know if to advance without resorting to equals geometry::equal_to comparator; - if (comparator(intersection_point, current)) + // Move backward and forward + RobustPoint direction_points[2]; + for (int i = 0; i < 2; i++) { - // It should be equal only once. But otherwise we skip it (in "add") - it = advance_circular(it, ring, robust_policy, current, seg_id, false); + int index = turn.operations[operation_index].index_in_robust_ring; + int piece_index = turn.operations[operation_index].piece_index; + while(comparator(pieces[piece_index].robust_ring[index], intersection_point)) + { + move_index(pieces, index, piece_index, i == 0 ? -1 : 1); + } + direction_points[i] = pieces[piece_index].robust_ring[index]; } - RobustPoint incoming = current; - - if (comparator(intersection_point, current)) - { - it = advance_circular(it, ring, robust_policy, current, real_seg_id); - } - else - { - // Don't upgrade the ID - it = advance_circular(it, ring, robust_policy, current, seg_id); - } - for (int defensive_check = 0; - comparator(intersection_point, current) && defensive_check < n; - defensive_check++) - { - it = advance_circular(it, ring, robust_policy, current, real_seg_id); - } - - info.add(incoming, current, intersection_point, turn_index, operation_index, real_seg_id); + info.add(direction_points[0], direction_points[1], intersection_point, + turn_index, operation_index, real_seg_id); } diff --git a/include/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp b/include/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp index cea2e2077..4f33bd40b 100644 --- a/include/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp +++ b/include/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp @@ -111,9 +111,11 @@ struct buffer_turn_operation : public detail::overlay::traversal_turn_operation { int piece_index; + int index_in_robust_ring; inline buffer_turn_operation() : piece_index(-1) + , index_in_robust_ring(-1) {} }; @@ -128,11 +130,13 @@ struct buffer_turn_info > { RobustPoint robust_point; + RobustPoint mapped_robust_point; // alas... we still need to adapt our points, offsetting them 1 integer to be co-located with neighbours bool is_opposite; intersection_location_type location; int count_within; + int count_on_offsetted; int count_on_occupied; int count_on_multi; @@ -144,6 +148,7 @@ struct buffer_turn_info : is_opposite(false) , location(location_ok) , count_within(0) + , count_on_offsetted(0) , count_on_occupied(0) , count_on_multi(0) {} diff --git a/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp b/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp index 92e2b90a5..1a0c4af4b 100644 --- a/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp +++ b/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp @@ -153,6 +153,7 @@ struct buffered_piece_collection struct robust_turn { int turn_index; + int operation_index; robust_point_type point; segment_identifier seg_id; segment_ratio_type fraction; @@ -185,6 +186,7 @@ struct buffered_piece_collection turn_vector_type m_turns; buffered_ring_collection > offsetted_rings; // indexed by multi_index + std::vector< std::vector < robust_point_type > > robust_offsetted_rings; buffered_ring_collection traversed_rings; segment_identifier current_segment_id; @@ -194,9 +196,6 @@ struct buffered_piece_collection { }; - typedef std::map > occupation_map_type; - occupation_map_type m_occupation_map; - struct redundant_turn { inline bool operator()(buffer_turn_info_type const& turn) const @@ -330,41 +329,75 @@ struct buffered_piece_collection int const geometry_code = detail::within::point_in_geometry(turn.robust_point, pc.robust_ring); - if (geometry_code == 1) + switch (geometry_code) { - turn.count_within++; + case 1 : turn.count_within++; break; + case 0 : turn.count_on_offsetted++; break; + } + + + } + + template + inline void adapt_mapped_robust_point(OccupationMap const& map, + buffer_turn_info_type& turn, int distance) const + { + for (int x = -distance; x <= distance; x++) + { + for (int y = -distance; y <= distance; y++) + { + robust_point_type rp = turn.robust_point; + geometry::set<0>(rp, geometry::get<0>(rp) + x); + geometry::set<1>(rp, geometry::get<1>(rp) + y); + if (map.find(rp) != map.end()) + { + turn.mapped_robust_point = rp; + return; + } + } } } - - // The "get_left_turn" process indicates, if it is a u/u turn (both only applicable - // for union, two separate turns), which is indicated in the map. If done so, set - // the other to "none", it is part of an occupied situation and should not be followed. - - inline void get_occupation() + inline void get_occupation(int distance = 0) { + typedef std::map + < + robust_point_type, + buffer_occupation_info, + geometry::less + > occupation_map_type; + + occupation_map_type occupation_map; + // 1: Add all intersection points to occupation map typedef typename boost::range_iterator::type const_iterator_type; typedef typename boost::range_iterator::type iterator_type; - for (const_iterator_type it = boost::begin(m_turns); + for (iterator_type it = boost::begin(m_turns); it != boost::end(m_turns); ++it) { - m_occupation_map[it->robust_point].m_count++; + if (it->count_on_offsetted >= 1) + { + if (distance > 0 && ! occupation_map.empty()) + { + adapt_mapped_robust_point(occupation_map, *it, distance); + } + occupation_map[it->mapped_robust_point].m_count++; + } } // 2: Remove all points from map which has only one - typename occupation_map_type::iterator it = m_occupation_map.begin(); - while (it != m_occupation_map.end()) + typename occupation_map_type::iterator it = occupation_map.begin(); + while (it != occupation_map.end()) { if (it->second.m_count <= 1) { typename occupation_map_type::iterator to_erase = it; ++it; - m_occupation_map.erase(to_erase); + occupation_map.erase(to_erase); } else { @@ -372,6 +405,11 @@ struct buffered_piece_collection } } + if (occupation_map.empty()) + { + return; + } + // 3: Add vectors (incoming->intersection-point, // intersection-point -> outgoing) // for all (co-located) points still present in the map @@ -390,19 +428,18 @@ struct buffered_piece_collection ++it, ++index) { typename occupation_map_type::iterator mit = - m_occupation_map.find(it->robust_point); + occupation_map.find(it->mapped_robust_point); - if (mit != m_occupation_map.end()) + if (mit != occupation_map.end()) { buffer_occupation_info& info = mit->second; // a: for (int i = 0; i < 2; i++) { - segment_identifier op_seg_id = it->operations[i].seg_id; - add_incoming_and_outgoing_angles(it->robust_point, - offsetted_rings[op_seg_id.multi_index], - m_robust_policy, - index, i, op_seg_id, info); + add_incoming_and_outgoing_angles(it->mapped_robust_point, *it, + m_pieces, + index, i, it->operations[i].seg_id, + info); } it->count_on_multi++; @@ -418,6 +455,24 @@ struct buffered_piece_collection } } + // X: Check rounding issues + if (distance == 0) + { + for (typename occupation_map_type::const_iterator it = occupation_map.begin(); + it != occupation_map.end(); ++it) + { + if (it->second.has_rounding_issues(it->first)) + { + std::cout << "Rounding issue! " << distance << std::endl; + if(distance == 0) + { + get_occupation(distance + 1); + return; + } + } + } + } + // 4: From these vectors, get the left turns // and mark all other turns as to-skip for (iterator_type it = boost::begin(m_turns); @@ -425,9 +480,9 @@ struct buffered_piece_collection ++it, ++index) { typename occupation_map_type::iterator mit = - m_occupation_map.find(it->robust_point); + occupation_map.find(it->mapped_robust_point); - if (mit != m_occupation_map.end()) + if (mit != occupation_map.end()) { buffer_occupation_info& info = mit->second; @@ -479,15 +534,12 @@ struct buffered_piece_collection for (typename boost::range_iterator::type it = boost::begin(m_turns); it != boost::end(m_turns); ++it) { - //if (it->count_on_occupied == 0) + typename std::vector::const_iterator pit; + for (pit = boost::begin(m_pieces); + pit != boost::end(m_pieces); + ++pit) { - typename std::vector::const_iterator pit; - for (pit = boost::begin(m_pieces); - pit != boost::end(m_pieces); - ++pit) - { - classify_turn(*it, *pit); - } + classify_turn(*it, *pit); } } } @@ -535,6 +587,30 @@ struct buffered_piece_collection } } + inline bool assert_indices_in_robust_rings() const + { + geometry::equal_to comparator; + for (typename boost::range_iterator::type it = + boost::begin(m_turns); it != boost::end(m_turns); ++it) + { + for (int i = 0; i < 2; i++) + { + robust_point_type const &p1 + = m_pieces[it->operations[i].piece_index].robust_ring + [it->operations[i].index_in_robust_ring]; + robust_point_type const &p2 = it->robust_point; + if (! comparator(p1, p2)) + { + std::cout << "FAILURE " << std::endl; + return false; + } + } + + + } + return true; + } + inline void rescale_pieces() { for (typename piece_vector_type::iterator it = boost::begin(m_pieces); @@ -584,11 +660,13 @@ struct buffered_piece_collection boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index) { geometry::recalculate(it->robust_point, it->point, m_robust_policy); + it->mapped_robust_point = it->robust_point; robust_turn turn; turn.turn_index = index; turn.point = it->robust_point; for (int i = 0; i < 2; i++) { + turn.operation_index = i; turn.seg_id = it->operations[i].seg_id; turn.fraction = it->operations[i].fraction; m_pieces[it->operations[i].piece_index].robust_turns.push_back(turn); @@ -613,22 +691,28 @@ struct buffered_piece_collection { std::sort(pc.robust_turns.begin(), pc.robust_turns.end(), buffer_operation_less()); } + // Walk through them, in reverse to insert at right index + int index_offset = pc.robust_turns.size() - 1; for (typename std::vector::const_reverse_iterator rit = pc.robust_turns.rbegin(); rit != pc.robust_turns.rend(); - ++rit) + ++rit, --index_offset) { - int const offsetted_count = pc.last_segment_index - pc.first_seg_id.segment_index; int const index_in_vector = 1 + rit->seg_id.segment_index - piece_segment_index; BOOST_ASSERT ( - index_in_vector > 0 && index_in_vector < 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++; + + m_turns[rit->turn_index].operations[rit->operation_index].index_in_robust_ring = index_in_vector + index_offset; } } } + + BOOST_ASSERT(assert_indices_in_robust_rings()); } template @@ -654,8 +738,9 @@ struct buffered_piece_collection rescale_pieces(); - get_occupation(); classify_turns(); + get_occupation(); + classify_inside(); check_remaining_points(input_geometry, distance_strategy); diff --git a/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp b/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp index 9183f8286..dcf9d1753 100644 --- a/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp +++ b/include/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp @@ -132,6 +132,7 @@ struct buffered_piece_collection_with_mapper out << " " << (it->count_within > 0 ? "w" : "") << (it->count_on_multi > 0 ? "m" : "") << (it->count_on_occupied > 0 ? "o" : "") + << (it->count_on_offsetted > 0 ? "b" : "") // offsetted border //<< it->debug_string ; @@ -141,13 +142,13 @@ struct buffered_piece_collection_with_mapper // out << it->operations[0].enriched.travels_to_vertex_index // << "/" << it->operations[1].enriched.travels_to_vertex_index; - offsets[it->robust_point] += 10; - int offset = offsets[it->robust_point]; + offsets[it->mapped_robust_point] += 10; + int offset = offsets[it->mapped_robust_point]; mapper.map(it->point, fill, 6); mapper.text(it->point, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:9px;", 5, offset); - offsets[it->robust_point] += 25; + offsets[it->mapped_robust_point] += 25; } } }