diff --git a/doc/release_notes.qbk b/doc/release_notes.qbk index 9717d430a..4e3ea1e99 100644 --- a/doc/release_notes.qbk +++ b/doc/release_notes.qbk @@ -20,6 +20,7 @@ [*Solved tickets] * [@https://svn.boost.org/trac/boost/ticket/9245 9245] Check for process errors in make_qbk.py +* [@https://svn.boost.org/trac/boost/ticket/9081 9081] Booleans create self-intersecting polygons from non-self-intersecting polygons [/=================] [heading Boost 1.55] diff --git a/extensions/test/algorithms/buffer/multi_point_buffer.cpp b/extensions/test/algorithms/buffer/multi_point_buffer.cpp index 25289ade6..fec6b810c 100644 --- a/extensions/test/algorithms/buffer/multi_point_buffer.cpp +++ b/extensions/test/algorithms/buffer/multi_point_buffer.cpp @@ -201,8 +201,8 @@ int test_main(int, char* []) //test_all >(); test_all >(); -#ifdef NDEBUG - // only in release mode + +#ifdef BOOST_GEOMETRY_BUFFER_TEST_GROWTH for (int i = 5; i <= 50; i++) { test_growth >(i, 20); 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 6da40c446..cbbc259ec 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp @@ -16,16 +16,14 @@ #include #include -#include -#include #include // Silence warning C4127: conditional expression is constant #if defined(_MSC_VER) -#pragma warning(push) -#pragma warning(disable : 4127) +#pragma warning(push) +#pragma warning(disable : 4127) #endif @@ -39,7 +37,7 @@ class turn_info_exception : public geometry::exception public: // NOTE: "char" will be replaced by enum in future version - inline turn_info_exception(char const method) + inline turn_info_exception(char const method) { message = "Boost.Geometry Turn exception: "; message += method; @@ -288,8 +286,14 @@ struct touch : public base_turn_handler int const side_pk_p = side.pk_wrt_p1(); int const side_qk_q = side.qk_wrt_q1(); + bool const both_continue = side_pk_p == 0 && side_qk_q == 0; + bool const robustness_issue_in_continue = both_continue && side_pk_q2 != 0; + bool const q_turns_left = side_qk_q == 1; - bool const block_q = side_qk_p1 == 0 && ! same(side_qi_p1, side_qk_q); + bool const block_q = side_qk_p1 == 0 + && ! same(side_qi_p1, side_qk_q) + && ! robustness_issue_in_continue + ; // If Pk at same side as Qi/Qk // (the "or" is for collinear case) @@ -479,20 +483,6 @@ struct equal : public base_turn_handler // oppositely if (side_pk_q2 == 0 && side_pk_p == side_qk_p) { - // After fixing robustness with integer: if they both continue collinearly, - // we might miss next intersection point because other rounding factor... - // Therefore we also look at the real-side (FP side) - typedef typename geometry::coordinate_type::type coordinate_type; - if (boost::is_floating_point::value) - { - coordinate_type sv_pk_p = strategy::side::side_by_triangle<>::side_value(pi, pj, pk); - coordinate_type sv_qk_p = strategy::side::side_by_triangle<>::side_value(pi, pj, qk); - if (sv_pk_p != sv_qk_p) - { - ui_else_iu(sv_pk_p > sv_qk_p, ti); - return; - } - } both(ti, operation_continue); return; @@ -587,6 +577,12 @@ struct collinear : public base_turn_handler - if Q arrives and Q turns left: union for Q (=intersection for P) - if Q arrives and Q turns right: intersection for Q (=union for P) + ROBUSTNESS: p and q are collinear, so you would expect + that side qk//p1 == pk//q1. But that is not always the case + in near-epsilon ranges. Then decision logic is different. + If p arrives, q is further, so the angle qk//p1 is (normally) + more precise than pk//p1 + */ template < @@ -620,6 +616,9 @@ struct collinear : public base_turn_handler : side_q ; + int const side_pk = side.pk_wrt_q1(); + int const side_qk = side.qk_wrt_p1(); + // See comments above, // resulting in a strange sort of mathematic rule here: // The arrival-info multiplied by the relevant side @@ -627,22 +626,15 @@ struct collinear : public base_turn_handler int const product = arrival * side_p_or_q; - if (product == 0 && (side_p != 0 || side_q != 0)) - { - // If q is collinear, p turns left - if (side_p != 0 && side_q == 0) - { - ui_else_iu(side_p == 1, ti); - return; - } - if (side_q != 0 && side_p == 0) - { - ui_else_iu(side_q == -1, ti); - return; - } - } + // 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(product == 0) + if (robustness_issue) + { + handle_robustness(ti, arrival, side_p, side_q, side_pk, side_qk); + } + else if(product == 0) { both(ti, operation_continue); } @@ -652,6 +644,37 @@ struct collinear : public base_turn_handler } } + static inline void handle_robustness(TurnInfo& ti, int arrival, + int side_p, int side_q, int side_pk, int side_qk) + { + // 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)) + { + use_p_for_union = false; + } + if (arrival == 1 && (consistent_side_p == 1 || consistent_side_q == -1)) + { + use_p_for_union = true; + } + + //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 @@ -693,9 +716,22 @@ private : typename IntersectionInfo > static inline bool set_tp(Point const& ri, Point const& rj, Point const& rk, int side_rk_r, - Point const& si, Point const& sj, + bool const handle_robustness, Point const& si, Point const& sj, int side_rk_s, TurnInfo& tp, IntersectionInfo const& intersection_info) { + if (handle_robustness) + { + // For Robustness: also calculate rk w.r.t. the other line. Because they are collinear opposite, that direction should be the reverse of the first direction. + // If this is not the case, we make it all-collinear, so zero + if (side_rk_r != 0 && side_rk_r != -side_rk_s) + { +#ifdef BOOST_GEOMETRY_DEBUG_ROBUSTNESS + std::cout << "Robustness correction: " << side_rk_r << " / " << side_rk_s << std::endl; +#endif + side_rk_r = 0; + } + } + operation_type blocked = operation_blocked; switch(side_rk_r) { @@ -764,7 +800,7 @@ public: // If P arrives within Q, there is a turn dependent on P if (dir_info.arrival[0] == 1 - && set_tp<0>(pi, pj, pk, side.pk_wrt_p1(), qi, qj, tp, intersection_info)) + && set_tp<0>(pi, pj, pk, side.pk_wrt_p1(), true, qi, qj, side.pk_wrt_q1(), tp, intersection_info)) { AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); *out++ = tp; @@ -772,7 +808,7 @@ public: // If Q arrives within P, there is a turn dependent on Q if (dir_info.arrival[1] == 1 - && set_tp<1>(qi, qj, qk, side.qk_wrt_q1(), pi, pj, tp, intersection_info)) + && set_tp<1>(qi, qj, qk, side.qk_wrt_q1(), false, pi, pj, side.qk_wrt_p1(), tp, intersection_info)) { AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); *out++ = tp; @@ -861,7 +897,7 @@ struct assign_null_policy static bool const include_degenerate = false; static bool const include_opposite = false; - template + template < typename Info, typename Point1, @@ -920,35 +956,12 @@ struct get_turn_info { typedef model::referring_segment segment_type1; typedef model::referring_segment segment_type2; - - typedef typename select_coordinate_type::type coordinate_type; - - typedef model::point - < - typename geometry::robust_type::type, - geometry::dimension::value, - typename geometry::coordinate_system::type - > robust_point_type; - - robust_point_type pi_rob, pj_rob, pk_rob, qi_rob, qj_rob, qk_rob; - geometry::zoom_to_robust(pi, pj, pk, qi, qj, qk, pi_rob, pj_rob, pk_rob, qi_rob, qj_rob, qk_rob); - - typedef geometry::strategy::side::side_by_triangle<> side; - side_info robust_sides; - robust_sides.set<0>(side::apply(qi_rob, qj_rob, pi_rob), - side::apply(qi_rob, qj_rob, pj_rob)); - robust_sides.set<1>(side::apply(pi_rob, pj_rob, qi_rob), - side::apply(pi_rob, pj_rob, qj_rob)); - - bool const p_equals = detail::equals::equals_point_point(pi_rob, pj_rob); - bool const q_equals = detail::equals::equals_point_point(qi_rob, qj_rob); - segment_type1 p1(pi, pj), p2(pj, pk); segment_type2 q1(qi, qj), q2(qj, qk); - side_calculator side_calc(pi_rob, pj_rob, pk_rob, qi_rob, qj_rob, qk_rob); + side_calculator side_calc(pi, pj, pk, qi, qj, qk); - typename strategy::return_type result = strategy::apply(p1, q1, robust_sides, p_equals, q_equals); + typename strategy::return_type result = strategy::apply(p1, q1); char const method = result.template get<1>().how; @@ -991,7 +1004,7 @@ struct get_turn_info else { // Swap p/q - side_calculator swapped_side_calc(qi_rob, qj_rob, qk_rob, pi_rob, pj_rob, pk_rob); + side_calculator swapped_side_calc(qi, qj, qk, pi, pj, pk); policy::template apply<1>(qi, qj, qk, pi, pj, pk, tp, result.template get<0>(), result.template get<1>(), swapped_side_calc); @@ -1110,7 +1123,7 @@ struct get_turn_info #if defined(_MSC_VER) -#pragma warning(pop) +#pragma warning(pop) #endif #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_GET_TURN_INFO_HPP diff --git a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp index 4c991d778..1da85cf2b 100644 --- a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/include/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -29,6 +29,10 @@ #include +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) +# include +#endif + namespace boost { namespace geometry { @@ -70,18 +74,6 @@ inline typename geometry::point_type::type get_from_index( } #endif -/*** -template -inline std::string rdebug(T const& value) -{ - if (math::equals(value, 0)) return "'0'"; - if (math::equals(value, 1)) return "'1'"; - if (value < 0) return "<0"; - if (value > 1) return ">1"; - return "<0..1>"; -} -***/ - /*! \see http://mathworld.wolfram.com/Line-LineIntersection.html */ @@ -105,7 +97,7 @@ struct relate_cartesian_segments static inline void debug_segments(std::string const& header, segment_type1 const& a, segment_type2 const& b) { std::cout << "Robustness issue: " << header << std::endl; - std::cout + std::cout << "A: " << wkt(a) << std::endl << "B: " << wkt(b) << std::endl ; @@ -129,13 +121,26 @@ struct relate_cartesian_segments coordinate_type const& dx_b, coordinate_type const& dy_b) { typedef side::side_by_triangle side; + side_info sides; coordinate_type const zero = 0; bool const a_is_point = math::equals(dx_a, zero) && math::equals(dy_a, zero); bool const b_is_point = math::equals(dx_b, zero) && math::equals(dy_b, zero); - // Get sides of a relative to b, and b relative to a - side_info sides; + if(a_is_point && b_is_point) + { + if(math::equals(get<1,0>(a), get<1,0>(b)) && math::equals(get<1,1>(a), get<1,1>(b))) + { + Policy::degenerate(a, true); + } + else + { + return Policy::disjoint(); + } + } + + 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) @@ -154,42 +159,19 @@ struct relate_cartesian_segments , detail::get_from_index<1>(a) , detail::get_from_index<1>(b)) ); - return apply(a, b, dx_a, dy_a, dx_b, dy_b, sides, a_is_point, b_is_point); - } - - static inline return_type apply(segment_type1 const& a, segment_type2 const& b, - side_info& sides, bool a_is_point, bool b_is_point) - { - coordinate_type const dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir - coordinate_type const dx_b = get<1, 0>(b) - get<0, 0>(b); - coordinate_type const dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir - coordinate_type const dy_b = get<1, 1>(b) - get<0, 1>(b); - return apply(a, b, dx_a, dy_a, dx_b, dy_b, sides, a_is_point, b_is_point); - } - - static inline return_type apply(segment_type1 const& a, segment_type2 const& b, - coordinate_type const& dx_a, coordinate_type const& dy_a, - coordinate_type const& dx_b, coordinate_type const& dy_b, - side_info& sides, bool a_is_point, bool b_is_point) - { - if(a_is_point && b_is_point) - { - if(math::equals(get<1,0>(a), get<1,0>(b)) && math::equals(get<1,1>(a), get<1,1>(b))) - { - Policy::degenerate(a, true); - } - else - { - return Policy::disjoint(); - } - } bool collinear = sides.collinear(); + robustness_verify_collinear(a, b, a_is_point, b_is_point, 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, are not disjoint @@ -235,13 +217,18 @@ struct relate_cartesian_segments { return Policy::disjoint(); } + + if (robustness_verify_disjoint_at_one_collinear(a, b, sides)) + { + return Policy::disjoint(); + } + } } if(collinear) { - // Use dimension 0 for collinear cases if differences in x exceeds differences in y - if (math::abs(dx_a) + math::abs(dx_b) >= math::abs(dy_a) + math::abs(dy_b)) + if (collinear_use_first) { return relate_collinear<0>(a, b); } @@ -317,6 +304,312 @@ private : return true; } + 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) + ; + } + + static inline void robustness_verify_collinear( + segment_type1 const& a, segment_type2 const& b, + bool a_is_point, bool b_is_point, + side_info& sides, + bool& collinear) + { + bool only_0_collinear = sides.zero<0>() && ! b_is_point && ! sides.zero<1>(); + bool only_1_collinear = sides.zero<1>() && ! a_is_point && ! sides.zero<0>(); + if (only_0_collinear || only_1_collinear) + { + typename geometry::point_type::type a0 = detail::get_from_index<0>(a); + typename geometry::point_type::type a1 = detail::get_from_index<1>(a); + typename geometry::point_type::type b0 = detail::get_from_index<0>(b); + typename geometry::point_type::type b1 = detail::get_from_index<1>(b); + bool ae = false, be = false; + + // If one of the segments is collinear, the other is probably so too. + side_info check; + coordinate_type factor = 1; + int iteration = 0; + bool collinear_consistent = false; + do + { + typedef side::side_by_triangle side; + + // We have a robustness issue. We keep increasing epsilon until we have a consistent set + coordinate_type const two = 2; + factor *= two; + coordinate_type epsilon = math::relaxed_epsilon(factor); + check.set<0> + ( + side::apply_with_epsilon(b0, b1, a0, epsilon), + side::apply_with_epsilon(b0, b1, a1, epsilon) + ); + check.set<1> + ( + side::apply_with_epsilon(a0, a1, b0, epsilon), + side::apply_with_epsilon(a0, a1, b1, epsilon) + ); + ae = point_equals_with_epsilon(a0, a1, epsilon); + be = point_equals_with_epsilon(b0, b1, epsilon); + + collinear_consistent = true; + if ( (check.zero<0>() && ! be && ! check.zero<1>()) + || (check.zero<1>() && ! ae && ! check.zero<0>()) + ) + { + collinear_consistent = false; + } + +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout + << "*** collinear_consistent: " + << iteration << std::boolalpha + << " consistent: " << collinear_consistent + << " equals: " << ae << "," << be + << " epsilon: " << epsilon + << " "; + check.debug(); +#endif + + + } while (! collinear_consistent && iteration++ < 10); + + sides = check; + collinear = sides.collinear(); + } + } + + 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) + ) + ) + { + + typename geometry::point_type::type a0 = detail::get_from_index<0>(a); + typename geometry::point_type::type a1 = detail::get_from_index<1>(a); + typename geometry::point_type::type b0 = detail::get_from_index<0>(b); + typename geometry::point_type::type b1 = detail::get_from_index<1>(b); + + side_info check; + coordinate_type factor = 1; + coordinate_type epsilon = math::relaxed_epsilon(factor); + int iteration = 1; + bool points_meet = false; + bool meeting_consistent = false; + do + { + typedef side::side_by_triangle side; + + // We have a robustness issue. We keep increasing epsilon until we have a consistent set + coordinate_type const two = 2; + factor *= two; + epsilon = math::relaxed_epsilon(factor); + check.set<0> + ( + side::apply_with_epsilon(b0, b1, a0, epsilon), + side::apply_with_epsilon(b0, b1, a1, epsilon) + ); + check.set<1> + ( + side::apply_with_epsilon(a0, a1, b0, epsilon), + side::apply_with_epsilon(a0, a1, b1, epsilon) + ); + + meeting_consistent = true; + if (check.meeting()) + { + points_meet = point_equals_with_epsilon + ( + select(check.zero_index<0>(), a), + select(check.zero_index<1>(), b), + epsilon + ); + if (! points_meet) + { + meeting_consistent = false; + + } + } + +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout + << "*** meeting_consistent: " + << iteration << std::boolalpha + << " consistent: " << meeting_consistent + << " epsilon: " << epsilon + << " "; + check.debug(); +#endif + + + } while (! meeting_consistent && iteration++ < 10); + + + sides = check; + + if (! sides.meeting() + && ((sides.zero<0>() && !sides.zero<1>()) + || (! sides.zero<0>() && sides.zero<1>()) + ) + ) + { + // Set sides to zero + sides.set<0>(0,0); + sides.set<1>(0,0); +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout << "ADAPTED New side info: " << std::endl; + sides.debug(); +#endif + } + + collinear = sides.collinear(); + + if (collinear_use_first && analyse_equal<0>(a, b)) + { + collinear_use_first = false; +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout << "Use [1] to check collinearity" << std::endl; +#endif + } + else if (! collinear_use_first && analyse_equal<1>(a, b)) + { + collinear_use_first = true; +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout << "Use [0] to check collinearity" << std::endl; +#endif + } + } + } + } + + // 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; + } + + 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 bool point_equals_with_epsilon(Point1 const& point1, Point2 const& point2, T const& epsilon) + { + // Check if difference is within espilon range (epsilon can be 0 for integer) + return math::abs(geometry::get<0>(point1) - geometry::get<0>(point2)) <= epsilon + && math::abs(geometry::get<1>(point1) - geometry::get<1>(point2)) <= epsilon + ; + } + + + // 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 verify_disjoint(segment_type1 const& a, segment_type2 const& b) @@ -336,8 +629,6 @@ private : 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); - - // TODO: these should be passed as integer too -> normal comparisons possible if (math::smaller(a_2, b_1) || math::larger(a_1, b_2)) //if (a_2 < b_1 || a_1 > b_2) { diff --git a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp index 9561a8ada..ed2cedcb6 100644 --- a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp +++ b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp @@ -65,7 +65,7 @@ public : return geometry::detail::determinant ( - dx, dy, + dx, dy, dpx, dpy ); @@ -99,8 +99,41 @@ public : promoted_type const s = side_value(p1, p2, p); promoted_type const zero = promoted_type(); - return math::equals(s, zero) ? 0 - : s > zero ? 1 + return math::equals(s, zero) ? 0 + : s > zero ? 1 + : -1; + } + + template + static inline int apply_with_epsilon(P1 const& p1, P2 const& p2, P const& p, T const& epsilon) + { + typedef typename boost::mpl::if_c + < + boost::is_void::type::value, + typename select_most_precise + < + typename select_most_precise + < + typename coordinate_type::type, + typename coordinate_type::type + >::type, + typename coordinate_type

::type + >::type, + CalculationType + >::type coordinate_type; + + // Promote float->double, small int->int + typedef typename select_most_precise + < + coordinate_type, + double + >::type promoted_type; + + promoted_type const s = side_value(p1, p2, p); + promoted_type const zero = promoted_type(); + + return math::abs(s - zero) < epsilon ? 0 + : s > zero ? 1 : -1; } }; diff --git a/include/boost/geometry/strategies/side_info.hpp b/include/boost/geometry/strategies/side_info.hpp index 3c2c1798a..d113eaa8f 100644 --- a/include/boost/geometry/strategies/side_info.hpp +++ b/include/boost/geometry/strategies/side_info.hpp @@ -16,8 +16,9 @@ #include #include -#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION -#include + +#if defined(BOOST_GEOMETRY_DEBUG_INTERSECTION) || defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) +# include #endif namespace boost { namespace geometry @@ -25,8 +26,8 @@ namespace boost { namespace geometry // Silence warning C4127: conditional expression is constant #if defined(_MSC_VER) -#pragma warning(push) -#pragma warning(disable : 4127) +#pragma warning(push) +#pragma warning(disable : 4127) #endif /*! @@ -145,13 +146,13 @@ public : return sides[Which].first == 0 ? 0 : 1; } -#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION +#if defined(BOOST_GEOMETRY_DEBUG_INTERSECTION) || defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) inline void debug() const { std::cout << sides[0].first << " " << sides[0].second << " " << sides[1].first << " " - << sides[1].second + << sides[1].second << std::endl; } #endif @@ -167,7 +168,7 @@ public : }; #if defined(_MSC_VER) -#pragma warning(pop) +#pragma warning(pop) #endif }} // namespace boost::geometry diff --git a/test/algorithms/difference.cpp b/test/algorithms/difference.cpp index a7b4bf533..df6d0c2a6 100644 --- a/test/algorithms/difference.cpp +++ b/test/algorithms/difference.cpp @@ -258,24 +258,19 @@ void test_all() 1, 0, 13); ***/ -#ifdef _MSC_VER -#ifdef TEST_ISOVIST test_one("isovist", isovist1[0], isovist1[1], - if_typed_tt(4, 2), 0, 0.279121891701124, - if_typed_tt(4, 3), 0, if_typed_tt(224.889211358929, 223.777), - if_typed_tt(0.001, 0.2)); + if_typed_tt(4, 2), -1, 0.279121, + 4, -1, 224.8892, + if_typed_tt(0.001, 0.1)); // SQL Server gives: 0.279121891701124 and 224.889211358929 // PostGIS gives: 0.279121991127244 and 224.889205853156 -#endif - test_one("ggl_list_20110306_javier", ggl_list_20110306_javier[0], ggl_list_20110306_javier[1], 1, -1, 71495.3331, 2, -1, 8960.49049); -#endif test_one("ggl_list_20110307_javier", ggl_list_20110307_javier[0], ggl_list_20110307_javier[1], diff --git a/test/algorithms/intersection.cpp b/test/algorithms/intersection.cpp index dbb1efc22..cedf52a7a 100644 --- a/test/algorithms/intersection.cpp +++ b/test/algorithms/intersection.cpp @@ -15,8 +15,6 @@ #include #include -#define TEST_ISOVIST - #include #include @@ -161,21 +159,14 @@ void test_areal() bool const open = bg::closure::value == bg::open; -#ifdef TEST_ISOVIST -#ifdef _MSC_VER test_one("isovist", isovist1[0], isovist1[1], - 1, if_typed(18, 20), 88.19203, + 1, 19, 88.19203, if_typed_tt(0.01, 0.1)); // SQL Server gives: 88.1920416352664 // PostGIS gives: 88.19203677911 -#endif -#endif - - //std::cout << typeid(ct).name() << std::endl; - if (! ccw && open) { // Pointcount for ttmath/double (both 5) or float (4) diff --git a/test/algorithms/overlay/get_turn_info.cpp b/test/algorithms/overlay/get_turn_info.cpp index 77eb77830..a750d1a12 100644 --- a/test/algorithms/overlay/get_turn_info.cpp +++ b/test/algorithms/overlay/get_turn_info.cpp @@ -323,14 +323,15 @@ void test_all() method_collinear, 5, 5, "ui"); // The next two cases are changed (BSG 2013-09-24), they contain turn info (#buffer_rt_g) + // In new approach they are changed back (BSG 2013-10-20) test_both("ccx1", 5, 1, 5, 6, 5, 8, // p 5, 5, 5, 7, 3, 8, // q - method_collinear, 5, 6, "iu"); // "cc"); + method_collinear, 5, 6, "cc"); // "iu"); test_both("cxc1", 5, 1, 5, 6, 7, 8, // p 5, 3, 5, 5, 5, 7, // q - method_collinear, 5, 5, "iu"); // "cc"); + method_collinear, 5, 5, "cc"); // "iu"); // Bug in case #54 of "overlay_cases.hpp" test_both("c_bug1", diff --git a/test/algorithms/overlay/robustness/ticket_9081.cpp b/test/algorithms/overlay/robustness/ticket_9081.cpp index 1e1aa97d5..7ce2e4e68 100644 --- a/test/algorithms/overlay/robustness/ticket_9081.cpp +++ b/test/algorithms/overlay/robustness/ticket_9081.cpp @@ -9,10 +9,7 @@ // Adapted from: the attachment of ticket 9081 -#define BOOST_GEOMETRY_DEBUG_HAS_SELF_INTERSECTIONS -#define BOOST_GEOMETRY_DEBUG_SEGMENT_IDENTIFIER -//#define BOOST_GEOMETRY_DEBUG_INTERSECTION -//#define CHECK_SELF_INTERSECTIONS +#define CHECK_SELF_INTERSECTIONS #include #include @@ -78,14 +75,16 @@ inline void debug_with_svg(int index, char method, Geometry const& a, Geometry c int main() { - int num_orig=50; - int num_rounds=20000; - //int num_rounds=16000; + int num_orig = 50; + int num_rounds = 20000; srand(1234); - std::cout< genesis; int pj; + + std::string wkt1, wkt2, operation; + try { @@ -116,13 +115,22 @@ int main() for(int j=0;j0) + if(boost::geometry::area(mp_i) > 0) { std::ostringstream out; out << "intersection(" << genesis[a] << " , " << genesis[b] << ")"; genesis[poly_list.size()] = out.str(); poly_list.push_back(mp_i); } - if(boost::geometry::area(mp_d)>0) + if(boost::geometry::area(mp_d) > 0) { std::ostringstream out; out << "difference(" << genesis[a] << " - " << genesis[b] << ")"; genesis[poly_list.size()] = out.str(); poly_list.push_back(mp_d); } - if(boost::geometry::area(mp_e)>0) + if(boost::geometry::area(mp_e) > 0) { std::ostringstream out; out << "difference(" << genesis[b] << ", " << genesis[a] << ")"; @@ -200,7 +208,11 @@ int main() } catch(std::exception const& e) { - std::cout << e.what() << " at " << pj << std::endl; + std::cout << e.what() + << " in " << operation << " at " << pj << std::endl + << wkt1 << std::endl + << wkt2 << std::endl + << std::endl; } catch(...) { diff --git a/test/algorithms/test_union.hpp b/test/algorithms/test_union.hpp index b8ccbca71..710b92101 100644 --- a/test/algorithms/test_union.hpp +++ b/test/algorithms/test_union.hpp @@ -1,4 +1,4 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) +// Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. @@ -37,16 +37,21 @@ template void test_union(std::string const& caseid, G1 const& g1, G2 const& g2, - std::size_t expected_count, int expected_hole_count, + int expected_count, int expected_hole_count, int expected_point_count, double expected_area, double percentage) { typedef typename bg::coordinate_type::type coordinate_type; std::vector clip; + +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout << "*** UNION " << caseid << std::endl; +#endif + bg::union_(g1, g2, clip); typename bg::default_area_result::type area = 0; - int n = 0; + std::size_t n = 0; std::size_t holes = 0; for (typename std::vector::iterator it = clip.begin(); it != clip.end(); ++it) @@ -56,6 +61,7 @@ void test_union(std::string const& caseid, G1 const& g1, G2 const& g2, n += bg::num_points(*it, true); } +#if ! defined(BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE) { // Test inserter functionality // Test if inserter returns output-iterator (using Boost.Range copy) @@ -78,35 +84,34 @@ void test_union(std::string const& caseid, G1 const& g1, G2 const& g2, BOOST_CHECK_EQUAL(boost::size(clip), boost::size(inserted) - 1); BOOST_CHECK_CLOSE(area_inserted, expected_area, percentage); } +#endif - /*** - std::cout << "case: " << caseid - << " n: " << n +#if defined(BOOST_GEOMETRY_DEBUG_ROBUSTNESS) + std::cout << "*** case: " << caseid << " area: " << area + << " points: " << n << " polygons: " << boost::size(clip) << " holes: " << holes << std::endl; - ***/ +#endif - if (expected_point_count >= 0) - { - BOOST_CHECK_MESSAGE(bg::math::abs(n - expected_point_count) < 3, - "union: " << caseid - << " #points expected: " << expected_point_count - << " detected: " << n - << " type: " << (type_for_assert_message()) - ); - } + BOOST_CHECK_MESSAGE(expected_point_count < 0 || std::abs(int(n) - expected_point_count) < 3, + "union: " << caseid + << " #points expected: " << expected_point_count + << " detected: " << n + << " type: " << (type_for_assert_message()) + ); - BOOST_CHECK_MESSAGE(clip.size() == expected_count, + BOOST_CHECK_MESSAGE(expected_count < 0 || int(clip.size()) == expected_count, "union: " << caseid << " #clips expected: " << expected_count << " detected: " << clip.size() << " type: " << (type_for_assert_message()) ); - BOOST_CHECK_MESSAGE(expected_hole_count < 0 || holes == std::size_t(expected_hole_count), + + BOOST_CHECK_MESSAGE(expected_hole_count < 0 || int(holes) == expected_hole_count, "union: " << caseid << " #holes expected: " << expected_hole_count << " detected: " << holes @@ -161,7 +166,7 @@ void test_union(std::string const& caseid, G1 const& g1, G2 const& g2, template void test_one(std::string const& caseid, std::string const& wkt1, std::string const& wkt2, - std::size_t expected_count, std::size_t expected_hole_count, + int expected_count, int expected_hole_count, int expected_point_count, double expected_area, double percentage = 0.001) { diff --git a/test/algorithms/union.cpp b/test/algorithms/union.cpp index 59b6be43f..42f5e6f61 100644 --- a/test/algorithms/union.cpp +++ b/test/algorithms/union.cpp @@ -15,8 +15,6 @@ #include #include -#define TEST_ISOVIST - #include #include @@ -216,48 +214,43 @@ void test_areal() test_one("ggl_list_20110306_javier", ggl_list_20110306_javier[0], ggl_list_20110306_javier[1], 1, 1, 16, 80456.4904910401); - + test_one("ggl_list_20110307_javier", ggl_list_20110307_javier[0], ggl_list_20110307_javier[1], 1, 1, 13, 20016.4); test_one("ggl_list_20110627_phillip", ggl_list_20110627_phillip[0], ggl_list_20110627_phillip[1], - 1, 0, - if_typed(5, if_typed_tt(8, 8)), + 1, 0, + if_typed(5, if_typed_tt(8, 8)), 14729.07145); - + // FP might return different amount of points test_one("ggl_list_20110716_enrico", ggl_list_20110716_enrico[0], ggl_list_20110716_enrico[1], 1, 1, - if_typed(18, if_typed(-1, 17)), + if_typed(18, if_typed(-1, 17)), 129904.197692871); - test_one("ggl_list_20110820_christophe", + test_one("ggl_list_20110820_christophe", ggl_list_20110820_christophe[0], ggl_list_20110820_christophe[1], - 2, // Previously it was 1 for double. The intersections are missed now, they fall in eps margins, arbitrary to generate them... - 0, - if_typed_tt(9, 8), + -1, // Either 1 or 2, depending if the intersection/turn point (eps.region) is missed + 0, + if_typed_tt(9, 8), 67.3550722317627); -#ifdef TEST_ISOVIST -#ifdef _MSC_VER test_one("isovist", isovist1[0], isovist1[1], 1, 0, - if_typed(72, if_typed(67, 73)), - 313.36036462, 0.1); + -1, + 313.36036462, 0.01); // SQL Server gives: 313.360374193241 // PostGIS gives: 313.360364623393 -#endif -#endif - // Ticket 5103 https://svn.boost.org/trac/boost/ticket/5103 // This ticket was actually reported for Boost.Polygon // We check it for Boost.Geometry as well. @@ -301,7 +294,7 @@ void test_areal() if (test_rt_i_rev) { test_one("buffer_rt_i_rev", buffer_rt_i[1], buffer_rt_i[0], - 1, 0, 13, 13.6569); + 1, 0, 13, 13.6569); } test_one("buffer_rt_j", buffer_rt_j[0], buffer_rt_j[1],