From d03d2e17a40919724899c7836a802f13748fec00 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Fri, 16 Mar 2012 16:58:26 +0000 Subject: [PATCH] [geometry] robustness fixes (all found by buffer robustness tests) [SVN r77351] --- .../strategies/cartesian/cart_intersect.hpp | 129 +++++++++++++----- .../boost/geometry/strategies/side_info.hpp | 20 +++ 2 files changed, 112 insertions(+), 37 deletions(-) diff --git a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp index 14f71d529..e20daaf3d 100644 --- a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -152,6 +152,33 @@ struct relate_cartesian_segments 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; + } + } + if (sides.same<0>() || sides.same<1>()) { // Both points are at same side of other segment, we can leave @@ -206,51 +233,47 @@ struct relate_cartesian_segments if (r < zero || r > one) { - if (sides.crossing() || sides.touching()) + 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) { - // ROBUSTNESS: the r value can in epsilon-cases be 1.14, while (with perfect arithmetic) - // it should be one. 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) - // TODO: find more cases - if (r > one) - { - // std::cout << "ROBUSTNESS: correction of r " << r << std::endl; - r = one; - } - else if (r < zero) - { - // std::cout << "ROBUSTNESS: correction of r " << r << std::endl; - r = zero; - } + r = one; + } + else if (r < zero) + { + r = zero; } - - if (r < zero) - { - if (r < -epsilon) - { - return Policy::disjoint(); - } - r = zero; - } - else if (r > one) - { - if (r > one + epsilon) - { - return Policy::disjoint(); - } - r = one; - } } } } if(collinear) { - if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_a)) + if (math::abs(dx_a) + math::abs(dx_b) < math::abs(dy_a) + math::abs(dy_b)) { // Y direction contains larger segments (maybe dx is zero) return relate_collinear<1>(a, b); @@ -268,6 +291,38 @@ struct relate_cartesian_segments private : + template + static inline bool double_check_disjoint(segment_type1 const& a, + segment_type2 const& b) + { + coordinate_type a_1, a_2, b_1, b_2; + bool a_swapped = false, b_swapped = false; + detail::segment_arrange(a, a_1, a_2, a_swapped); + detail::segment_arrange(b, b_1, b_2, b_swapped); + return math::smaller(a_2, b_1) || math::larger(a_1, b_2); + } + + template + static inline typename point_type::type select(int index, Segment const& segment) + { + return index == 0 + ? detail::get_from_index<0>(segment) + : detail::get_from_index<1>(segment) + ; + } + + // 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_equals(Point1 const& point1, Point2 const& point2) + { + return math::equals(get<0>(point1), get<0>(point2)) + && math::equals(get<1>(point1), get<1>(point2)) + ; + } + + + template static inline return_type relate_collinear(segment_type1 const& a, segment_type2 const& b) diff --git a/include/boost/geometry/strategies/side_info.hpp b/include/boost/geometry/strategies/side_info.hpp index a28c1914c..8744e2e40 100644 --- a/include/boost/geometry/strategies/side_info.hpp +++ b/include/boost/geometry/strategies/side_info.hpp @@ -81,12 +81,32 @@ public : && sides[1].second == 0 && sides[0].second == 0); } + inline bool meeting() const + { + // Two of them (in each segment) zero, two not + return one_zero<0>() && one_zero<1>(); + } + template inline bool zero() const { return sides[Which].first == 0 && sides[Which].second == 0; } + template + inline bool one_zero() const + { + return (sides[Which].first == 0 && sides[Which].second != 0) + || (sides[Which].first != 0 && sides[Which].second == 0); + } + + template + inline int zero_index() const + { + return sides[Which].first == 0 ? 0 : 1; + } + + inline void debug() const { std::cout << sides[0].first << " "