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 e90910c9bc..e6fb7cd782 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp @@ -15,7 +15,6 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_GET_TURN_INFO_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_GET_TURN_INFO_HPP - #include #include @@ -28,9 +27,6 @@ #include #include -#include - -#include #include // Silence warning C4127: conditional expression is constant @@ -131,15 +127,7 @@ struct base_turn_handler return 0; } - typedef typename select_coordinate_type - < - typename UniqueSubRange1::point_type, - typename UniqueSubRange2::point_type - >::type coordinate_type; - - typedef detail::distance_measure dm_type; - - dm_type const dm = get_distance_measure(range_p.at(range_index), range_p.at(range_index + 1), range_q.at(point_index)); + auto const dm = get_distance_measure(range_p.at(range_index), range_p.at(range_index + 1), range_q.at(point_index)); return dm.measure == 0 ? 0 : dm.measure > 0 ? 1 : -1; } @@ -186,7 +174,8 @@ struct base_turn_handler { ti.method = method; - // For touch/touch interior always take the intersection point 0 (there is only one). + // For touch/touch interior always take the intersection point 0 + // (usually there is only one - but if collinear is handled as touch, both could be taken). static int const index = 0; geometry::convert(info.intersections[index], ti.point); @@ -225,9 +214,8 @@ struct base_turn_handler { // TODO: use comparable distance for point-point instead - but that // causes currently cycling include problems - typedef typename geometry::coordinate_type::type ctype; - ctype const dx = get<0>(a) - get<0>(b); - ctype const dy = get<1>(a) - get<1>(b); + auto const dx = get<0>(a) - get<0>(b); + auto const dy = get<1>(a) - get<1>(b); return dx * dx + dy * dy; } @@ -259,19 +247,10 @@ struct base_turn_handler // pk/q2 is considered as collinear, but there might be // a tiny measurable difference. If so, use that. // Calculate pk // qj-qk - typedef detail::distance_measure - < - typename select_coordinate_type - < - typename UniqueSubRange1::point_type, - typename UniqueSubRange2::point_type - >::type - > dm_type; - - const bool p_closer = + bool const p_closer = ti.operations[IndexP].remaining_distance < ti.operations[IndexQ].remaining_distance; - dm_type const dm + auto const dm = p_closer ? get_distance_measure(range_q.at(index_q - 1), range_q.at(index_q), range_p.at(index_p)) @@ -282,8 +261,7 @@ struct base_turn_handler { // Not truely collinear, distinguish for union/intersection // If p goes left (positive), take that for a union - - bool p_left = p_closer ? dm.is_positive() : dm.is_negative(); + bool const p_left = p_closer ? dm.is_positive() : dm.is_negative(); ti.operations[IndexP].operation = p_left ? operation_union : operation_intersection; @@ -347,14 +325,9 @@ struct touch_interior : public base_turn_handler // Therefore handle it as a normal touch (two segments arrive at the // intersection point). It currently checks for zero, but even a // distance a little bit larger would do. - typedef typename geometry::coordinate_type - < - typename UniqueSubRange::point_type - >::type coor_t; - - coor_t const location = distance_measure(info.intersections[0], non_touching_range.at(1)); - coor_t const zero = 0; - bool const result = math::equals(location, zero); + auto const dm = distance_measure(info.intersections[0], non_touching_range.at(1)); + decltype(dm) const zero = 0; + bool const result = math::equals(dm, zero); return result; } @@ -540,16 +513,8 @@ struct touch : public base_turn_handler // >----->P qj is LEFT of P1 and pi is LEFT of Q2 // (the other way round is also possible) - typedef typename select_coordinate_type - < - typename UniqueSubRange1::point_type, - typename UniqueSubRange2::point_type - >::type coordinate_type; - - typedef detail::distance_measure dm_type; - - dm_type const dm_qj_p1 = get_distance_measure(range_p.at(0), range_p.at(1), range_q.at(1)); - dm_type const dm_pi_q2 = get_distance_measure(range_q.at(1), range_q.at(2), range_p.at(0)); + auto const dm_qj_p1 = get_distance_measure(range_p.at(0), range_p.at(1), range_q.at(1)); + auto const dm_pi_q2 = get_distance_measure(range_q.at(1), range_q.at(2), range_p.at(0)); if (dm_qj_p1.measure > 0 && dm_pi_q2.measure > 0) { @@ -564,8 +529,8 @@ struct touch : public base_turn_handler return true; } - dm_type const dm_pj_q1 = get_distance_measure(range_q.at(0), range_q.at(1), range_p.at(1)); - dm_type const dm_qi_p2 = get_distance_measure(range_p.at(1), range_p.at(2), range_q.at(0)); + auto const dm_pj_q1 = get_distance_measure(range_q.at(0), range_q.at(1), range_p.at(1)); + auto const dm_qi_p2 = get_distance_measure(range_p.at(1), range_p.at(2), range_q.at(0)); if (dm_pj_q1.measure > 0 && dm_qi_p2.measure > 0) { @@ -813,17 +778,9 @@ struct equal : public base_turn_handler // They turn to the same side, or continue both collinearly // Without rescaling, to check for union/intersection, // try to check side values (without any thresholds) - typedef typename select_coordinate_type - < - typename UniqueSubRange1::point_type, - typename UniqueSubRange2::point_type - >::type coordinate_type; - - typedef detail::distance_measure dm_type; - - dm_type const dm_pk_q2 + auto const dm_pk_q2 = get_distance_measure(range_q.at(1), range_q.at(2), range_p.at(2)); - dm_type const dm_qk_p2 + auto const dm_qk_p2 = get_distance_measure(range_p.at(1), range_p.at(2), range_q.at(2)); if (dm_qk_p2.measure != dm_pk_q2.measure) @@ -965,8 +922,57 @@ template > struct collinear : public base_turn_handler { + template + < + typename IntersectionInfo, + typename UniqueSubRange1, + typename UniqueSubRange2, + typename DirInfo + > + static bool handle_as_equal(IntersectionInfo const& info, + UniqueSubRange1 const& range_p, + UniqueSubRange2 const& range_q, + DirInfo const& dir_info) + { +#if defined(BOOST_GEOMETRY_USE_RESCALING) + return false; +#endif + int const arrival_p = dir_info.arrival[0]; + int const arrival_q = dir_info.arrival[1]; + if (arrival_p * arrival_q != -1 || info.count != 2) + { + // Code below assumes that either p or q arrives in the other segment + return false; + } + + auto const dm = distance_measure(info.intersections[1], + arrival_p == 1 ? range_q.at(1) : range_p.at(1)); + decltype(dm) const zero = 0; + return math::equals(dm, zero); + } + /* - arrival P pk//p1 qk//q1 product* case result + Either P arrives within Q (arrival_p == -1) or Q arrives within P. + + Typical situation: + ^q2 + / + ^q1 + / ____ ip[1] + //|p1 } this section of p/q is colllinear + q0// | } ____ ip[0] + / | + / v + p0 p2 + + P arrives (at p1) in segment Q (between q0 and q1). + Therefore arrival_p == 1 + P (p2) goes to the right (-1). Follow P for intersection, or follow Q for union. + Therefore if (arrival_p==1) and side_p==-1, result = iu + + Complete table: + + arrival P pk//p1 qk//q1 product case result 1 1 1 CLL1 ui -1 1 -1 CLL2 iu 1 1 1 CLR1 ui @@ -980,7 +986,9 @@ struct collinear : public base_turn_handler 1 0 0 CC1 cc -1 0 0 CC2 cc - *product = arrival * (pk//p1 or qk//q1) + Resulting in the rule: + The arrival-info multiplied by the relevant side delivers the result. + product = arrival * (pk//p1 or qk//q1) Stated otherwise: - if P arrives: look at turn P @@ -989,13 +997,6 @@ struct collinear : public base_turn_handler - if P arrives and P turns right: intersection for P - 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 < @@ -1016,9 +1017,9 @@ struct collinear : public base_turn_handler // Copy the intersection point in TO direction assign_point(ti, method_collinear, info, non_opposite_to_index(info)); - int const arrival = dir_info.arrival[0]; + int const arrival_p = dir_info.arrival[0]; // Should not be 0, this is checked before - BOOST_GEOMETRY_ASSERT(arrival != 0); + BOOST_GEOMETRY_ASSERT(arrival_p != 0); bool const has_pk = ! range_p.is_last_segment(); bool const has_qk = ! range_q.is_last_segment(); @@ -1026,19 +1027,15 @@ struct collinear : public base_turn_handler int const side_q = has_qk ? side.qk_wrt_q1() : 0; // If p arrives, use p, else use q - int const side_p_or_q = arrival == 1 + int const side_p_or_q = arrival_p == 1 ? side_p : side_q ; - // See comments above, - // resulting in a strange sort of mathematic rule here: - // The arrival-info multiplied by the relevant side - // delivers a consistent result. + // Calculate product according to comments above. + int const product = arrival_p * side_p_or_q; - int const product = arrival * side_p_or_q; - - if(product == 0) + if (product == 0) { both(ti, operation_continue); } @@ -1186,11 +1183,11 @@ private : { TurnInfo tp = tp_model; - int const p_arrival = info.d_info().arrival[0]; - int const q_arrival = info.d_info().arrival[1]; + int const arrival_p = info.d_info().arrival[0]; + int const arrival_q = info.d_info().arrival[1]; // If P arrives within Q, there is a turn dependent on P - if ( p_arrival == 1 + if ( arrival_p == 1 && ! range_p.is_last_segment() && set_tp<0>(side.pk_wrt_p1(), tp, info.i_info()) ) { @@ -1200,7 +1197,7 @@ private : } // If Q arrives within P, there is a turn dependent on Q - if ( q_arrival == 1 + if ( arrival_q == 1 && ! range_q.is_last_segment() && set_tp<1>(side.qk_wrt_q1(), tp, info.i_info()) ) { @@ -1212,8 +1209,8 @@ private : if (AssignPolicy::include_opposite) { // Handle cases not yet handled above - if ((q_arrival == -1 && p_arrival == 0) - || (p_arrival == -1 && q_arrival == 0)) + if ((arrival_q == -1 && arrival_p == 0) + || (arrival_p == -1 && arrival_q == 0)) { for (unsigned int i = 0; i < 2; i++) { @@ -1420,9 +1417,10 @@ struct get_turn_info // Collinear if ( ! inters.d_info().opposite ) { - - if ( inters.d_info().arrival[0] == 0 ) + if (inters.d_info().arrival[0] == 0 + || collinear::handle_as_equal(inters.i_info(), range_p, range_q, inters.d_info())) { + // Both segments arrive at the second intersection point handle_as_equal = true; } else diff --git a/test/algorithms/buffer/buffer_multi_polygon.cpp b/test/algorithms/buffer/buffer_multi_polygon.cpp index f5421ec44b..94535088b8 100644 --- a/test/algorithms/buffer/buffer_multi_polygon.cpp +++ b/test/algorithms/buffer/buffer_multi_polygon.cpp @@ -349,6 +349,34 @@ static std::string const nores_wt_1 static std::string const nores_wt_2 = "MULTIPOLYGON(((1 1,2 2,2 1,1 1)),((0 2,0 3,1 3,0 2)),((4 3,5 4,5 3,4 3)),((4 3,3 3,4 4,4 3)))"; +static std::string const nores_b8e6 + = "MULTIPOLYGON(((4 4,5 5,5 4,4 4)),((4 2,5 2,5 1,3 1,4 2)),((3 1,4 0,3 0,3 1)))"; +static std::string const nores_2881 + = "MULTIPOLYGON(((2 5,2 6,3 5,2 5)),((1 7,0 7,0 8,1 8,1 7)),((1 7,1 6,0 6,1 7)))"; +static std::string const nores_3af0 + = "MULTIPOLYGON(((1 8,0 8,0 9,1 9,1 8)),((1 8,1 7,0 7,1 8)),((2 4,3 5,3 4,2 4)),((2 6,2 7,3 6,2 6)))"; +static std::string const nores_6061 + = "MULTIPOLYGON(((2 8,2 9,3 8,2 8)),((4 3,4 4,5 4,4 3)),((7 2,7 3,8 2,7 2)),((5 3,6 4,6 3,5 3)),((2 6,3 7,3 6,2 6)),((2 3,3 2,3 1,2 1,2 3)),((2 3,3 4,3 3,2 3)))"; + +static std::string const nores_1ea1 + = "MULTIPOLYGON(((2 0,2 1,3 0,2 0)),((7 5,6 4,5 3,5 4,5 5,7 5)),((2 3,1 3,1 4,2 3)),((2 3,2 4,3 3,2 3)))"; +static std::string const nores_804e + = "MULTIPOLYGON(((4 8,4 9,5 8,4 8)),((3 9,3 10,4 10,3 9)),((0 7,0 8,1 7,0 7)),((4 6,3 6,3 7,4 6)),((4 6,4 7,5 7,4 6)))"; +static std::string const nores_37f6 + = "MULTIPOLYGON(((4 1,5 2,5 1,4 1)),((1 0,1 1,2 1,2 0,1 0)),((0 3,0 4,1 4,1 3,0 3)),((2 2,2 3,3 2,2 2)))"; + +static std::string const nores_495d + = "MULTIPOLYGON(((2 0,2 1,3 0,2 0)),((2 3,3 4,3 3,2 3)),((5 1,5 2,6 2,5 1)),((4 3,4 2,3 2,4 3)))"; + +// rescaled +static std::string const res_ebc4 + = "MULTIPOLYGON(((3 9,3 10,4 9,3 9)),((9 5,9 6,10 6,10 5,9 5)),((8 8,8 9,9 9,8 8)),((4 8,3 7,3 8,4 8)),((4 8,5 9,6 9,6 8,4 8)),((4 5,3 4,3 5,4 5)),((4 5,5 6,5 5,4 5)))"; +static std::string const res_8618 + = "MULTIPOLYGON(((6 2,7 3,7 2,6 2)),((4 3,5 4,5 3,4 3)),((3 0,3 1,4 0,3 0)),((8 7,8 8,9 8,8 7)),((0 7,0 8,1 8,1 7,0 7)),((2 2,1 2,1 3,1 4,2 4,2 2)),((2 2,2 1,1 1,2 2)))"; +static std::string const res_3b4d + = "MULTIPOLYGON(((8 0,9 1,9 0,8 0)),((3 4,2 4,2 5,2 6,3 6,3 4)),((3 4,4 3,3 3,3 4)),((3 8,3 9,4 9,3 8)),((0 5,0 6,1 6,0 5)),((7 3,8 4,8 3,7 3)),((5 5,6 6,6 5,5 5)))"; + + static std::string const neighbouring = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)),((10 10,10 20,20 20,20 10,10 10)))"; @@ -483,7 +511,6 @@ void test_all() test_one("rt_p3", rt_p3, join_miter, end_flat, 22.3995, 1.0); test_one("rt_p4", rt_p4, join_miter, end_flat, 33.0563, 1.0); test_one("rt_p5", rt_p5, join_miter, end_flat, 17.0, 1.0); - test_one("rt_p6", rt_p6, join_miter, end_flat, 18.4853, 1.0); test_one("rt_p7", rt_p7, join_miter, end_flat, 26.2279, 1.0); test_one("rt_p8", rt_p8, join_miter, end_flat, 29.0563, 1.0); @@ -571,6 +598,26 @@ void test_all() test_one("nores_wt_1", nores_wt_1, join_round32, end_flat, 80.1609, 1.0); test_one("nores_wt_2", nores_wt_2, join_round32, end_flat, 22.1102, 1.0); + test_one("nores_b8e6", nores_b8e6, join_round32, end_flat, 19.8528, 1.0); + + test_one("nores_2881", nores_2881, join_round32, end_flat, 16.5517, 1.0); + test_one("nores_6061", nores_6061, join_round32, end_flat, 39.7371, 1.0); + test_one("nores_37f6", nores_37f6, join_round32, end_flat, 26.5339, 1.0); + +#if defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES) + // Cases not yet solved without rescaling (out of 241) + test_one("nores_1ea1", nores_1ea1, join_round32, end_flat, 28.9755, 1.0); + test_one("nores_804e", nores_804e, join_round32, end_flat, 26.4503, 1.0); + test_one("nores_3af0", nores_3af0, join_round32, end_flat, 22.1008, 1.0); + test_one("nores_495d", nores_495d, join_round32, end_flat, 23.4376, 1.0); +#endif + +#if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES) + // Erroneous cases with rescaling (out of 8) + test_one("res_ebc4", res_ebc4, join_round32, end_flat, 43.8877, 1.0); + test_one("res_8618", res_8618, join_round32, end_flat, 48.1085, 1.0); + test_one("res_3b4d", res_3b4d, join_round32, end_flat, 48.4739, 1.0); +#endif test_one("neighbouring_small", neighbouring, @@ -615,7 +662,7 @@ int test_main(int, char* []) #endif #if defined(BOOST_GEOMETRY_TEST_FAILURES) - BoostGeometryWriteExpectedFailures(1, 1, 2, 3); + BoostGeometryWriteExpectedFailures(4, 6, 3, 8); #endif return 0; diff --git a/test/algorithms/set_operations/difference/difference.cpp b/test/algorithms/set_operations/difference/difference.cpp index ef6b2e140f..6d3d49999d 100644 --- a/test/algorithms/set_operations/difference/difference.cpp +++ b/test/algorithms/set_operations/difference/difference.cpp @@ -624,7 +624,7 @@ int test_main(int, char* []) #if defined(BOOST_GEOMETRY_TEST_FAILURES) // Not yet fully tested for float and long double. // The difference algorithm can generate (additional) slivers - BoostGeometryWriteExpectedFailures(10, 11, 24, 16); + BoostGeometryWriteExpectedFailures(10, 11, 24, 15); #endif return 0;