From 2aed65855f76afd5f4beaf3429168efc895305c8 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Sun, 15 Apr 2012 11:44:15 +0000 Subject: [PATCH] [geometry] fix of several robustness issues in cart_intersect and get_turn_info found by testing buffer algorithm. Also restructured cart_intersect such that all robustness issues are handled in separate methods (could be policy later). Finally fixed ever circling iterator with range (for assignment) [SVN r77988] --- .../detail/overlay/get_turn_info.hpp | 73 ++-- .../iterators/ever_circling_iterator.hpp | 134 +++++-- .../strategies/cartesian/cart_intersect.hpp | 375 +++++++++++++----- .../boost/geometry/strategies/side_info.hpp | 32 ++ 4 files changed, 432 insertions(+), 182 deletions(-) diff --git a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp index 609825e34..b5725eda2 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp @@ -571,6 +571,9 @@ struct collinear : public base_turn_handler : side_q ; + int const side_pk = SideStrategy::apply(qi, qj, pk); + int const side_qk = SideStrategy::apply(pi, pj, qk); + // See comments above, // resulting in a strange sort of mathematic rule here: // The arrival-info multiplied by the relevant side @@ -578,53 +581,55 @@ struct collinear : public base_turn_handler int const product = arrival * side_p_or_q; - if(product == 0) + // Robustness: side_p is supposed to be equal to side_pk (because p/q are collinear) + // and side_q to side_qk + bool const robustness_issue = side_pk != side_p || side_qk != side_q; + + if (robustness_issue) + { + handle_robustness(ti, arrival, side_p, side_q, side_pk, side_qk); + } + else if(product == 0) { both(ti, operation_continue); } else { - int const side_pk = SideStrategy::apply(qi, qj, pk); - int const side_qk = SideStrategy::apply(pi, pj, qk); - - if (side_pk != side_p || side_qk != side_q) - { - //std::cout << "ROBUSTNESS -> Collinear " - // << " arr: " << arrival - // << " prod: " << product - // << " dir: " << side_p << " " << side_q - // << " rev: " << side_pk << " " << side_qk - // << std::endl; - - handle_robustness(ti, arrival, product, - side_p, side_q, side_pk, side_qk); - } - else - { - // The normal case - ui_else_iu(product == 1, ti); - } + ui_else_iu(product == 1, ti); } } - static inline void handle_robustness(TurnInfo& ti, - int arrival, int product, - int side_p, int side_q, - int side_pk, int side_qk) + static inline void handle_robustness(TurnInfo& ti, int arrival, + int side_p, int side_q, int side_pk, int side_qk) { - bool take_ui = product == 1; - if (product == arrival) + // We take the longer one, i.e. if q arrives in p (arrival == -1), + // then p exceeds q and we should take p for a union... + + bool use_p_for_union = arrival == -1; + + // ... unless one of the sides consistently directs to the other side + int const consistent_side_p = side_p == side_pk ? side_p : 0; + int const consistent_side_q = side_q == side_qk ? side_q : 0; + if (arrival == -1 && (consistent_side_p == -1 || consistent_side_q == 1)) { - if ((product == 1 && side_p == 1 && side_pk != 1) - || (product == -1 && side_q == 1 && side_qk != 1)) - { - //std::cout << "ROBUSTNESS: -> Reverse" << std::endl; - take_ui = ! take_ui; - } + use_p_for_union = false; + } + if (arrival == 1 && (consistent_side_p == 1 || consistent_side_q == -1)) + { + use_p_for_union = true; } - ui_else_iu(take_ui, ti); + //std::cout << "ROBUSTNESS -> Collinear " + // << " arr: " << arrival + // << " dir: " << side_p << " " << side_q + // << " rev: " << side_pk << " " << side_qk + // << " cst: " << cside_p << " " << cside_q + // << std::boolalpha << " " << use_p_for_union + // << std::endl; + + ui_else_iu(use_p_for_union, ti); } + }; template diff --git a/include/boost/geometry/iterators/ever_circling_iterator.hpp b/include/boost/geometry/iterators/ever_circling_iterator.hpp index 03921096e..566669e26 100644 --- a/include/boost/geometry/iterators/ever_circling_iterator.hpp +++ b/include/boost/geometry/iterators/ever_circling_iterator.hpp @@ -95,67 +95,115 @@ private: bool m_skip_first; }; - - template -class ever_circling_range_iterator - : public boost::iterator_adaptor - < - ever_circling_range_iterator, - typename boost::range_iterator::type - > +struct ever_circling_range_iterator + : public boost::iterator_facade + < + ever_circling_range_iterator, + typename boost::range_value::type const, + boost::random_access_traversal_tag + > { -public : - typedef typename boost::range_iterator::type iterator_type; + /// Constructor including the range it is based on + explicit inline ever_circling_range_iterator(Range& range) + : m_range(&range) + , m_iterator(boost::begin(range)) + , m_size(boost::size(range)) + , m_index(0) + {} - explicit inline ever_circling_range_iterator(Range& range, - bool skip_first = false) - : m_range(range) - , m_skip_first(skip_first) + /// Default constructor + explicit inline ever_circling_range_iterator() + : m_range(NULL) + , m_size(0) + , m_index(0) + {} + + inline ever_circling_range_iterator& operator=(ever_circling_range_iterator const& source) { - this->base_reference() = boost::begin(m_range); + m_range = source.m_range; + m_iterator = source.m_iterator; + m_size = source.m_size; + m_index = source.m_index; + return *this; } - explicit inline ever_circling_range_iterator(Range& range, iterator_type start, - bool skip_first = false) - : m_range(range) - , m_skip_first(skip_first) - { - this->base_reference() = start; - } - - /// Navigate to a certain position, should be in [start .. end], if at end - /// it will circle again. - inline void moveto(iterator_type it) - { - this->base_reference() = it; - check_end(); - } + typedef std::ptrdiff_t difference_type; private: - friend class boost::iterator_core_access; - inline void increment(bool possibly_skip = true) + inline typename boost::range_value::type const& dereference() const { - (this->base_reference())++; - check_end(possibly_skip); + return *m_iterator; } - inline void check_end(bool possibly_skip = true) + inline difference_type distance_to(ever_circling_range_iterator const& other) const { - if (this->base_reference() == boost::end(m_range)) + return other.m_index - this->m_index; + } + + inline bool equal(ever_circling_range_iterator const& other) const + { + return this->m_range == other.m_range + && this->m_index == other.m_index; + } + + inline void increment() + { + ++m_index; + if (m_index >= 0 && m_index < m_size) { - this->base_reference() = boost::begin(m_range); - if (m_skip_first && possibly_skip) - { - increment(false); - } + ++m_iterator; + } + else + { + update_iterator(); } } - Range& m_range; - bool m_skip_first; + inline void decrement() + { + --m_index; + if (m_index >= 0 && m_index < m_size) + { + --m_iterator; + } + else + { + update_iterator(); + } + } + + inline void advance(difference_type n) + { + if (m_index >= 0 && m_index < m_size + && m_index + n >= 0 && m_index + n < m_size) + { + m_index += n; + m_iterator += n; + } + else + { + m_index += n; + update_iterator(); + } + } + + inline void update_iterator() + { + while (m_index < 0) + { + m_index += m_size; + } + m_index = m_index % m_size; + this->m_iterator = boost::begin(*m_range) + m_index; + } + + Range* m_range; + typename boost::range_iterator::type m_iterator; + difference_type m_size; + difference_type m_index; }; diff --git a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp index e20daaf3d..ed0f50087 100644 --- a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -120,6 +120,8 @@ struct relate_cartesian_segments typedef side::side_by_triangle side; side_info sides; + bool collinear_use_first = math::abs(dx_a) + math::abs(dx_b) >= math::abs(dy_a) + math::abs(dy_b); + sides.set<0> ( side::apply(detail::get_from_index<0>(b) @@ -141,48 +143,16 @@ struct relate_cartesian_segments bool collinear = sides.collinear(); - if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>())) - { - // If one of the segments is collinear, the other must be as well. - // So handle it as collinear. - // (In float/double epsilon margins it can easily occur that one or two of them are -1/1) - // sides.debug(); - sides.set<0>(0,0); - sides.set<1>(0,0); - collinear = true; - } - - if (sides.meeting()) - { - // If two segments meet each other at their segment-points, two sides are zero, - // the other two are not (unless collinear but we don't mean those here). - // However, in near-epsilon ranges it can happen that two sides are zero - // but they do not meet at their segment-points. - // In that case they are nearly collinear and handled as such. - if (! point_equals - ( - select(sides.zero_index<0>(), a), - select(sides.zero_index<1>(), b) - ) - ) - { - //std::cout << "ROBUSTNESS: Suspicious "; - //sides.debug(); - //std::cout << std::endl; - //std::cout << " " << d << " " << da << std::endl; - //std::cout << " -> " << geometry::wkt(a) << " " << geometry::wkt(b) << std::endl; - //std::cout << " -> " << dx_a << " " << dy_a << " " << dx_b << " " << dy_b << std::endl; - - sides.set<0>(0,0); - sides.set<1>(0,0); - collinear = true; - } - } + robustness_verify_collinear(a, b, sides, collinear); + robustness_verify_meeting(a, b, sides, collinear, collinear_use_first); if (sides.same<0>() || sides.same<1>()) { // Both points are at same side of other segment, we can leave - return Policy::disjoint(); + if (robustness_verify_same_side(a, b, sides)) + { + return Policy::disjoint(); + } } // Degenerate cases: segments of single point, lying on other segment, non disjoint @@ -225,62 +195,31 @@ struct relate_cartesian_segments { r = da / d; - // Ratio should lie between 0 and 1 - // Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear - promoted_type const zero = 0; - promoted_type const one = 1; - promoted_type const epsilon = std::numeric_limits::epsilon(); - - if (r < zero || r > one) + if (! robustness_verify_r(a, b, r)) { - if (double_check_disjoint<0>(a, b) - || double_check_disjoint<1>(a, b)) - { - // Can still be disjoint (even if not one is left or right from another) - // This is e.g. in case #snake4 of buffer test. - return Policy::disjoint(); - } - - //std::cout << "ROBUSTNESS: correction of r " << r << std::endl; - // sides.debug() - - // ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic) - // it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting) - - // If segments are crossing (we can see that with the sides) - // and one is inside the other, there must be an intersection point. - // We correct for that. - // This is (only) in case #ggl_list_20110820_christophe in unit tests - - // If segments are touching (two sides zero), of course they should intersect - // This is (only) in case #buffer_rt_i in the unit tests) - - // If one touches in the middle, they also should intersect (#buffer_rt_j) - - // Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575 - - if (r > one) - { - r = one; - } - else if (r < zero) - { - r = zero; - } + return Policy::disjoint(); } + + robustness_handle_meeting(a, b, sides, dx_a, dy_a, wx, wy, d, r); + + if (robustness_verify_disjoint_at_one_collinear(a, b, sides)) + { + return Policy::disjoint(); + } + } } if(collinear) { - if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_b)) + if (collinear_use_first) { - // Y direction contains larger segments (maybe dx is zero) - return relate_collinear<1>(a, b); + return relate_collinear<0>(a, b); } else { - return relate_collinear<0>(a, b); + // Y direction contains larger segments (maybe dx is zero) + return relate_collinear<1>(a, b); } } @@ -291,8 +230,210 @@ struct relate_cartesian_segments private : + + // Ratio should lie between 0 and 1 + // Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear + template + static inline bool robustness_verify_r( + segment_type1 const& a, segment_type2 const& b, + T& r) + { + T const zero = 0; + T const one = 1; + if (r < zero || r > one) + { + if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) + { + // Can still be disjoint (even if not one is left or right from another) + // This is e.g. in case #snake4 of buffer test. + return false; + } + + //std::cout << "ROBUSTNESS: correction of r " << r << std::endl; + // sides.debug(); + + // ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic) + // it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting) + + // If segments are crossing (we can see that with the sides) + // and one is inside the other, there must be an intersection point. + // We correct for that. + // This is (only) in case #ggl_list_20110820_christophe in unit tests + + // If segments are touching (two sides zero), of course they should intersect + // This is (only) in case #buffer_rt_i in the unit tests) + + // If one touches in the middle, they also should intersect (#buffer_rt_j) + + // Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575 + + if (r > one) + { + r = one; + } + else if (r < zero) + { + r = zero; + } + } + return true; + } + + static inline void robustness_verify_collinear( + segment_type1 const& a, segment_type2 const& b, + side_info& sides, + bool& collinear) + { + if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>())) + { + // If one of the segments is collinear, the other must be as well. + // So handle it as collinear. + // (In float/double epsilon margins it can easily occur that one or two of them are -1/1) + // sides.debug(); + sides.set<0>(0,0); + sides.set<1>(0,0); + collinear = true; + } + } + + static inline void robustness_verify_meeting( + segment_type1 const& a, segment_type2 const& b, + side_info& sides, + bool& collinear, bool& collinear_use_first) + { + if (sides.meeting()) + { + // If two segments meet each other at their segment-points, two sides are zero, + // the other two are not (unless collinear but we don't mean those here). + // However, in near-epsilon ranges it can happen that two sides are zero + // but they do not meet at their segment-points. + // In that case they are nearly collinear and handled as such. + if (! point_equals + ( + select(sides.zero_index<0>(), a), + select(sides.zero_index<1>(), b) + ) + ) + { + sides.set<0>(0,0); + sides.set<1>(0,0); + collinear = true; + + if (collinear_use_first && analyse_equal<0>(a, b)) + { + collinear_use_first = false; + } + else if (! collinear_use_first && analyse_equal<1>(a, b)) + { + collinear_use_first = true; + } + + } + } + } + + // Verifies and if necessary correct missed touch because of robustness + // This is the case at multi_polygon_buffer unittest #rt_m + static inline bool robustness_verify_same_side( + segment_type1 const& a, segment_type2 const& b, + side_info& sides) + { + int corrected = 0; + if (sides.one_touching<0>()) + { + if (point_equals( + select(sides.zero_index<0>(), a), + select(0, b) + )) + { + sides.correct_to_zero<1, 0>(); + corrected = 1; + } + if (point_equals + ( + select(sides.zero_index<0>(), a), + select(1, b) + )) + { + sides.correct_to_zero<1, 1>(); + corrected = 2; + } + } + else if (sides.one_touching<1>()) + { + if (point_equals( + select(sides.zero_index<1>(), b), + select(0, a) + )) + { + sides.correct_to_zero<0, 0>(); + corrected = 3; + } + if (point_equals + ( + select(sides.zero_index<1>(), b), + select(1, a) + )) + { + sides.correct_to_zero<0, 1>(); + corrected = 4; + } + } + + return corrected == 0; + } + + static inline bool robustness_verify_disjoint_at_one_collinear( + segment_type1 const& a, segment_type2 const& b, + side_info const& sides) + { + if (sides.one_of_all_zero()) + { + if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) + { + return true; + } + } + return false; + } + + + // If r is one, or zero, segments should meet and their endpoints. + // Robustness issue: check if this is really the case. + // It turns out to be no problem, see buffer test #rt_s1 (and there are many cases generated) + // It generates an "ends in the middle" situation which is correct. + template + static inline void robustness_handle_meeting(segment_type1 const& a, segment_type2 const& b, + side_info& sides, + T const& dx_a, T const& dy_a, T const& wx, T const& wy, + T const& d, R const& r) + { + return; + + T const db = geometry::detail::determinant(dx_a, dy_a, wx, wy); + + R const zero = 0; + R const one = 1; + if (math::equals(r, zero) || math::equals(r, one)) + { + R rb = db / d; + if (rb <= 0 || rb >= 1 || math::equals(rb, 0) || math::equals(rb, 1)) + { + if (sides.one_zero<0>() && ! sides.one_zero<1>()) // or vice versa + { +#if defined(BOOST_GEOMETRY_COUNT_INTERSECTION_EQUAL) + extern int g_count_intersection_equal; + g_count_intersection_equal++; +#endif + sides.debug(); + std::cout << "E r=" << r << " r.b=" << rb << " "; + } + } + } + } + template - static inline bool double_check_disjoint(segment_type1 const& a, + static inline bool verify_disjoint(segment_type1 const& a, segment_type2 const& b) { coordinate_type a_1, a_2, b_1, b_2; @@ -321,7 +462,30 @@ private : ; } + // We cannot use geometry::equals here. Besides that this will be changed + // to compare segment-coordinate-values directly (not necessary to retrieve point first) + template + static inline bool point_equality(Point1 const& point1, Point2 const& point2, + bool& equals_0, bool& equals_1) + { + equals_0 = math::equals(get<0>(point1), get<0>(point2)); + equals_1 = math::equals(get<1>(point1), get<1>(point2)); + return equals_0 && equals_1; + } + template + static inline bool analyse_equal(segment_type1 const& a, segment_type2 const& b) + { + coordinate_type const a_1 = geometry::get<0, Dimension>(a); + coordinate_type const a_2 = geometry::get<1, Dimension>(a); + coordinate_type const b_1 = geometry::get<0, Dimension>(b); + coordinate_type const b_2 = geometry::get<1, Dimension>(b); + return math::equals(a_1, b_1) + || math::equals(a_2, b_1) + || math::equals(a_1, b_2) + || math::equals(a_2, b_2) + ; + } template static inline return_type relate_collinear(segment_type1 const& a, @@ -367,25 +531,28 @@ private : // Handle "equal", in polygon neighbourhood comparisons a common case - // Check if segments are equal... - bool const a1_eq_b1 = math::equals(get<0, 0>(a), get<0, 0>(b)) - && math::equals(get<0, 1>(a), get<0, 1>(b)); - bool const a2_eq_b2 = math::equals(get<1, 0>(a), get<1, 0>(b)) - && math::equals(get<1, 1>(a), get<1, 1>(b)); - if (a1_eq_b1 && a2_eq_b2) + bool const opposite = a_swapped ^ b_swapped; + bool const both_swapped = a_swapped && b_swapped; + + // Check if segments are equal or opposite equal... + bool const swapped_a1_eq_b1 = math::equals(a_1, b_1); + bool const swapped_a2_eq_b2 = math::equals(a_2, b_2); + + if (swapped_a1_eq_b1 && swapped_a2_eq_b2) { - return Policy::segment_equal(a, false); + return Policy::segment_equal(a, opposite); } - // ... or opposite equal - bool const a1_eq_b2 = math::equals(get<0, 0>(a), get<1, 0>(b)) - && math::equals(get<0, 1>(a), get<1, 1>(b)); - bool const a2_eq_b1 = math::equals(get<1, 0>(a), get<0, 0>(b)) - && math::equals(get<1, 1>(a), get<0, 1>(b)); - if (a1_eq_b2 && a2_eq_b1) - { - return Policy::segment_equal(a, true); - } + bool const swapped_a2_eq_b1 = math::equals(a_2, b_1); + bool const swapped_a1_eq_b2 = math::equals(a_1, b_2); + + bool const a1_eq_b1 = both_swapped ? swapped_a2_eq_b2 : a_swapped ? swapped_a2_eq_b1 : b_swapped ? swapped_a1_eq_b2 : swapped_a1_eq_b1; + bool const a2_eq_b2 = both_swapped ? swapped_a1_eq_b1 : a_swapped ? swapped_a1_eq_b2 : b_swapped ? swapped_a2_eq_b1 : swapped_a2_eq_b2; + + bool const a1_eq_b2 = both_swapped ? swapped_a2_eq_b1 : a_swapped ? swapped_a2_eq_b2 : b_swapped ? swapped_a1_eq_b1 : swapped_a1_eq_b2; + bool const a2_eq_b1 = both_swapped ? swapped_a1_eq_b2 : a_swapped ? swapped_a1_eq_b1 : b_swapped ? swapped_a2_eq_b2 : swapped_a2_eq_b1; + + // The rest below will return one or two intersections. @@ -393,7 +560,7 @@ private : // For IM it is important to know which relates to which. So this information is given, // without performance penalties to intersection calculation - bool const has_common_points = a1_eq_b1 || a1_eq_b2 || a2_eq_b1 || a2_eq_b2; + bool const has_common_points = swapped_a1_eq_b1 || swapped_a1_eq_b2 || swapped_a2_eq_b1 || swapped_a2_eq_b2; // "Touch" -> one intersection point -> one but not two common points @@ -413,7 +580,7 @@ private : // #4: a2<----a1 b1--->b2 (no arrival at all) // Where the arranged forms have two forms: // a_1-----a_2/b_1-------b_2 or reverse (B left of A) - if (has_common_points && (math::equals(a_2, b_1) || math::equals(b_2, a_1))) + if ((swapped_a2_eq_b1 || swapped_a1_eq_b2) && ! swapped_a1_eq_b1 && ! swapped_a2_eq_b2) { if (a2_eq_b1) return Policy::collinear_touch(get<1, 0>(a), get<1, 1>(a), 0, -1); if (a1_eq_b2) return Policy::collinear_touch(get<0, 0>(a), get<0, 1>(a), -1, 0); @@ -473,7 +640,6 @@ private : if (a1_eq_b1) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, arrival_a, -arrival_a, false); } - bool const opposite = a_swapped ^ b_swapped; // "Inside", a completely within b or b completely within a @@ -535,7 +701,6 @@ private : the picture might seem wrong but it (supposed to be) is right. */ - bool const both_swapped = a_swapped && b_swapped; if (b_1 < a_2 && a_2 < b_2) { // Left column, from bottom to top @@ -557,7 +722,7 @@ private : ; } // Nothing should goes through. If any we have made an error - // Robustness: it can occur here... + // std::cout << "Robustness issue, non-logical behaviour" << std::endl; return Policy::error("Robustness issue, non-logical behaviour"); } }; diff --git a/include/boost/geometry/strategies/side_info.hpp b/include/boost/geometry/strategies/side_info.hpp index 8744e2e40..f3a9da0df 100644 --- a/include/boost/geometry/strategies/side_info.hpp +++ b/include/boost/geometry/strategies/side_info.hpp @@ -44,6 +44,19 @@ public : sides[Which].second = second; } + template + inline void correct_to_zero() + { + if (Index == 0) + { + sides[Which].first = 0; + } + else + { + sides[Which].second = 0; + } + } + template inline int get() const { @@ -81,6 +94,15 @@ public : && sides[1].second == 0 && sides[0].second == 0); } + template + inline bool one_touching() const + { + // This is normally a situation which can't occur: + // If one is completely left or right, the other cannot touch + return one_zero() + && sides[1 - Which].first * sides[1 - Which].second == 1; + } + inline bool meeting() const { // Two of them (in each segment) zero, two not @@ -100,6 +122,16 @@ public : || (sides[Which].first != 0 && sides[Which].second == 0); } + inline bool one_of_all_zero() const + { + int const sum = std::abs(sides[0].first) + + std::abs(sides[0].second) + + std::abs(sides[1].first) + + std::abs(sides[1].second); + return sum == 3; + } + + template inline int zero_index() const {