diff --git a/include/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp b/include/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp index a35be052a..287a06b08 100644 --- a/include/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp @@ -441,18 +441,21 @@ inline void enrich_intersection_points(Turns& turns, if (turn.both(detail::overlay::operation_none) || turn.both(opposite_operation) + || turn.both(detail::overlay::operation_blocked) || (detail::overlay::is_self_turn(turn) && ! turn.is_clustered() && ! turn.both(target_operation))) { + // For all operations, discard xx and none/none // For intersections, remove uu to avoid the need to travel // a union (during intersection) in uu/cc clusters (e.g. #31,#32,#33) + // The ux is necessary to indicate impossible paths + // (especially if rescaling is removed) - // Similarly, for union, discard ii + // Similarly, for union, discard ii and ix - // Only keep self-uu-turns or self-ii-turns + // For self-turns, only keep uu / ii - // Blocked (or combination with blocked is still needed for difference) turn.discarded = true; turn.cluster_id = -1; continue; 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 edb05f97b..5cb5f78c7 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp @@ -85,7 +85,7 @@ struct base_turn_handler return side1 * side2 == 1; } - // Both continue + // Both get the same operation template static inline void both(TurnInfo& ti, operation_type const op) { @@ -131,6 +131,17 @@ struct base_turn_handler ? 1 : 0; } + template + static inline typename geometry::coordinate_type::type + distance_measure(Point1 const& a, Point2 const& b) + { + // 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); + return dx * dx + dy * dy; + } }; @@ -186,6 +197,9 @@ struct touch_interior : public base_turn_handler int const side_qk_q = has_k ? side.qk_wrt_q1() : 0; + // Only necessary if rescaling is turned off: + int const side_pj_q2 = side.pj_wrt_q2(); + if (side_qi_p == -1 && side_qk_p == -1 && side_qk_q == 1) { // Q turns left on the right side of P (test "MR3") @@ -195,10 +209,21 @@ struct touch_interior : public base_turn_handler } else if (side_qi_p == 1 && side_qk_p == 1 && side_qk_q == -1) { - // Q turns right on the left side of P (test "ML3") - // Union: take both operation - // Intersection: skip - both(ti, operation_union); + if (side_pj_q2 == -1) + { + // Q turns right on the left side of P (test "ML3") + // Union: take both operations + // Intersection: skip + both(ti, operation_union); + } + else + { + // q2 is collinear with p1, so it does not turn back. This + // can happen in floating point precision. In this case, + // block one of the operations to avoid taking that path. + ti.operations[index_p].operation = operation_union; + ti.operations[index_q].operation = operation_blocked; + } ti.touch_only = true; } else if (side_qi_p == side_qk_p && side_qi_p == side_qk_q) @@ -208,6 +233,29 @@ struct touch_interior : public base_turn_handler // Union: take left turn (Q if Q turns left, P if Q turns right) // Intersection: other turn unsigned int index = side_qk_q == 1 ? index_q : index_p; + if (side_pj_q2 == 0) + { + // Even though sides xk w.r.t. 1 are distinct, pj is collinear + // with q. Therefore swap the path + index = 1 - index; + } + + if (opposite(side_pj_q2, side_qi_p)) + { + // Without rescaling, floating point requires extra measures + int const side_qj_p1 = side.qj_wrt_p1(); + int const side_qj_p2 = side.qj_wrt_p2(); + + if (same(side_qj_p1, side_qj_p2)) + { + int const side_pj_q1 = side.pj_wrt_q1(); + if (opposite(side_pj_q1, side_pj_q2)) + { + index = 1 - index; + } + } + } + ti.operations[index].operation = operation_union; ti.operations[1 - index].operation = operation_intersection; ti.touch_only = true; @@ -220,10 +268,16 @@ struct touch_interior : public base_turn_handler // Collinearly in the same direction // (Q comes from left of P and turns left, // OR Q comes from right of P and turns right) - // Omit intersection point. + // Omit second intersection point. // Union: just continue // Intersection: just continue both(ti, operation_continue); + + // Calculate remaining distance. + // Q arrives at p, at point qj, so use qk for q + // and use pj for p + ti.operations[index_p].remaining_distance = distance_measure(ti.point, range_p.at(1)); + ti.operations[index_q].remaining_distance = distance_measure(ti.point, range_q.at(2)); } else { @@ -642,18 +696,6 @@ struct collinear : public base_turn_handler ? distance_measure(ti.point, range_q.at(2)) : distance_measure(ti.point, range_q.at(1)); } - - template - static inline typename geometry::coordinate_type::type - distance_measure(Point1 const& a, Point2 const& b) - { - // 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); - return dx * dx + dy * dy; - } }; template diff --git a/include/boost/geometry/algorithms/detail/overlay/get_turn_info_helpers.hpp b/include/boost/geometry/algorithms/detail/overlay/get_turn_info_helpers.hpp index c87c92e1d..087ca8060 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info_helpers.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info_helpers.hpp @@ -73,6 +73,12 @@ struct side_calculator inline int pk_wrt_q2() const { return m_side_strategy.apply(get_qj(), get_qk(), get_pk()); } inline int qk_wrt_p2() const { return m_side_strategy.apply(get_pj(), get_pk(), get_qk()); } + // Necessary when rescaling turns off: + inline int qj_wrt_p1() const { return m_side_strategy.apply(get_pi(), get_pj(), get_qj()); } + inline int qj_wrt_p2() const { return m_side_strategy.apply(get_pj(), get_pk(), get_qj()); } + inline int pj_wrt_q1() const { return m_side_strategy.apply(get_qi(), get_qj(), get_pj()); } + inline int pj_wrt_q2() const { return m_side_strategy.apply(get_qj(), get_qk(), get_pj()); } + inline point1_type const& get_pi() const { return m_range_p.at(0); } inline point1_type const& get_pj() const { return m_range_p.at(1); } inline point1_type const& get_pk() const { return m_range_p.at(2); } diff --git a/include/boost/geometry/algorithms/detail/overlay/traversal.hpp b/include/boost/geometry/algorithms/detail/overlay/traversal.hpp index 95f57941b..3a7d82ce0 100644 --- a/include/boost/geometry/algorithms/detail/overlay/traversal.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/traversal.hpp @@ -105,6 +105,23 @@ template > struct traversal { +private : + struct linked_turn_op_info + { + explicit linked_turn_op_info(signed_size_type ti = -1, int oi = -1, + signed_size_type nti = -1) + : turn_index(ti) + , op_index(oi) + , next_turn_index(nti) + , rank_index(-1) + {} + + signed_size_type turn_index; + int op_index; + signed_size_type next_turn_index; + signed_size_type rank_index; + }; + static const operation_type target_operation = operation_from_overlay::value; typedef typename side_compare::type side_compare_type; @@ -118,6 +135,7 @@ struct traversal point_type, SideStrategy, side_compare_type > sbs_type; +public : inline traversal(Geometry1 const& geometry1, Geometry2 const& geometry2, Turns& turns, Clusters const& clusters, RobustPolicy const& robust_policy, SideStrategy const& strategy, @@ -354,35 +372,72 @@ struct traversal // If both are valid candidates, take the one with minimal remaining // distance (important for #mysql_23023665 in buffer). - // Initialize with 0, automatically assigned on first result + signed_size_type next[2] = {0}; + bool possible[2] = {0}; + bool close[2] = {0}; + + for (int i = 0; i < 2; i++) + { + next[i] = turn.operations[i].enriched.get_next_turn_index(); + possible[i] = traverse_possible(next[i]); + close[i] = possible[i] && next[i] == start_turn_index; + } + + if (close[0] != close[1]) + { + // One of the operations will finish the ring. Take that one. + selected_op_index = close[0] ? 0 : 1; + debug_traverse(turn, turn.operations[selected_op_index], "Candidate cc closing"); + return true; + } + + if (OverlayType == overlay_buffer && possible[0] && possible[1]) + { + // Buffers sometimes have multiple overlapping pieces, where remaining + // distance could lead to the wrong choice. Take the matching operation. + + bool is_target[2] = {0}; + for (int i = 0; i < 2; i++) + { + turn_operation_type const& next_op = m_turns[next[i]].operations[i]; + is_target[i] = next_op.operation == target_operation; + } + + if (is_target[0] != is_target[1]) + { + // Take the matching operation + selected_op_index = is_target[0] ? 0 : 1; + debug_traverse(turn, turn.operations[selected_op_index], "Candidate cc target"); + return true; + } + } + + static bool const is_union = target_operation == operation_union; + typename turn_operation_type::comparable_distance_type - min_remaining_distance = 0; + best_remaining_distance = 0; bool result = false; for (int i = 0; i < 2; i++) { - turn_operation_type const& op = turn.operations[i]; - - signed_size_type const next_turn_index = op.enriched.get_next_turn_index(); - - if (! traverse_possible(next_turn_index)) + if (!possible[i]) { continue; } + turn_operation_type const& op = turn.operations[i]; + if (! result - || next_turn_index == start_turn_index - || op.remaining_distance < min_remaining_distance) + || (is_union && op.remaining_distance > best_remaining_distance) + || (!is_union && op.remaining_distance < best_remaining_distance)) { debug_traverse(turn, op, "First candidate cc", ! result); - debug_traverse(turn, op, "Candidate cc override (start)", - result && next_turn_index == start_turn_index); debug_traverse(turn, op, "Candidate cc override (remaining)", - result && op.remaining_distance < min_remaining_distance); + result && op.remaining_distance < best_remaining_distance); selected_op_index = i; - min_remaining_distance = op.remaining_distance; + best_remaining_distance = op.remaining_distance; result = true; } } @@ -712,24 +767,163 @@ struct traversal return false; } - inline bool select_turn_from_cluster(signed_size_type& turn_index, + inline signed_size_type get_rank(sbs_type const& sbs, + linked_turn_op_info const& info) const + { + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + typename sbs_type::rp const& rp = sbs.m_ranked_points[i]; + if (rp.turn_index == info.turn_index + && rp.operation_index == info.op_index + && rp.direction == sort_by_side::dir_to) + { + return rp.rank; + } + } + return -1; + } + + // Function checks simple cases, such as a cluster with two turns, + // arriving at the first turn, first turn points to second turn, + // second turn points further. + inline bool select_turn_from_cluster_linked(signed_size_type& turn_index, int& op_index, - signed_size_type start_turn_index, int start_op_index, + std::set const& ids, segment_identifier const& previous_seg_id) const { - bool const is_union = target_operation == operation_union; + typedef typename std::set::const_iterator sit_type; - turn_type const& turn = m_turns[turn_index]; - BOOST_ASSERT(turn.is_clustered()); + std::vector possibilities; + std::vector blocked; + for (sit_type it = ids.begin(); it != ids.end(); ++it) + { + signed_size_type cluster_turn_index = *it; + turn_type const& cluster_turn = m_turns[cluster_turn_index]; + if (cluster_turn.discarded) + { + continue; + } + if (is_self_turn(cluster_turn) + || cluster_turn.both(target_operation)) + { + // Not (yet) supported, can be cluster of u/u turns + return false; + } + for (int i = 0; i < 2; i++) + { + turn_operation_type const& op = cluster_turn.operations[i]; + turn_operation_type const& other_op = cluster_turn.operations[1 - i]; + signed_size_type const ni = op.enriched.get_next_turn_index(); + if (op.operation == target_operation + || op.operation == operation_continue) + { + if (ni == cluster_turn_index) + { + // Not (yet) supported, traveling to itself, can be + // hole + return false; + } + possibilities.push_back( + linked_turn_op_info(cluster_turn_index, i, ni)); + } + else if (op.operation == operation_blocked + && ! (ni == other_op.enriched.get_next_turn_index()) + && ids.count(ni) == 0) + { + // Points to turn, not part of this cluster, + // and that way is blocked. But if the other operation + // points at the same turn, it is still fine. + blocked.push_back( + linked_turn_op_info(cluster_turn_index, i, ni)); + } + } + } - typename Clusters::const_iterator mit = m_clusters.find(turn.cluster_id); - BOOST_ASSERT(mit != m_clusters.end()); + typedef typename std::vector::const_iterator const_it_type; - cluster_info const& cinfo = mit->second; - std::set const& ids = cinfo.turn_indices; + if (! blocked.empty()) + { + sbs_type sbs(m_strategy); - sbs_type sbs(m_strategy); + if (! fill_sbs(sbs, turn_index, ids, previous_seg_id)) + { + return false; + } + for (typename std::vector::iterator it = possibilities.begin(); + it != possibilities.end(); ++it) + { + linked_turn_op_info& info = *it; + info.rank_index = get_rank(sbs, info); + } + for (typename std::vector::iterator it = blocked.begin(); + it != blocked.end(); ++it) + { + linked_turn_op_info& info = *it; + info.rank_index = get_rank(sbs, info); + } + + + for (const_it_type it = possibilities.begin(); + it != possibilities.end(); ++it) + { + linked_turn_op_info const& lti = *it; + for (const_it_type bit = blocked.begin(); + bit != blocked.end(); ++bit) + { + linked_turn_op_info const& blti = *bit; + if (blti.next_turn_index == lti.next_turn_index + && blti.rank_index == lti.rank_index) + { + return false; + } + } + } + } + + // Traversal can either enter the cluster in the first turn, + // or it can start halfway. + // If there is one (and only one) possibility pointing outside + // the cluster, take that one. + linked_turn_op_info target; + for (const_it_type it = possibilities.begin(); + it != possibilities.end(); ++it) + { + linked_turn_op_info const& lti = *it; + if (ids.count(lti.next_turn_index) == 0) + { + if (target.turn_index >= 0 + && target.next_turn_index != lti.next_turn_index) + { + // Points to different target + return false; + } + if (OverlayType == overlay_buffer && target.turn_index > 0) + { + // Target already assigned, so there are more targets + // or more ways to the same target + return false; + } + + target = lti; + } + } + if (target.turn_index < 0) + { + return false; + } + + turn_index = target.turn_index; + op_index = target.op_index; + + return true; + } + + inline bool fill_sbs(sbs_type& sbs, + signed_size_type turn_index, + std::set const& ids, + segment_identifier const& previous_seg_id) const + { for (typename std::set::const_iterator sit = ids.begin(); sit != ids.end(); ++sit) { @@ -755,7 +949,39 @@ struct traversal { return false; } + turn_type const& turn = m_turns[turn_index]; sbs.apply(turn.point); + return true; + } + + + inline bool select_turn_from_cluster(signed_size_type& turn_index, + int& op_index, + signed_size_type start_turn_index, int start_op_index, + segment_identifier const& previous_seg_id) const + { + bool const is_union = target_operation == operation_union; + + turn_type const& turn = m_turns[turn_index]; + BOOST_ASSERT(turn.is_clustered()); + + typename Clusters::const_iterator mit = m_clusters.find(turn.cluster_id); + BOOST_ASSERT(mit != m_clusters.end()); + + cluster_info const& cinfo = mit->second; + std::set const& ids = cinfo.turn_indices; + + if (select_turn_from_cluster_linked(turn_index, op_index, ids, previous_seg_id)) + { + return true; + } + + sbs_type sbs(m_strategy); + + if (! fill_sbs(sbs, turn_index, ids, previous_seg_id)) + { + return false; + } bool result = false; diff --git a/include/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp b/include/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp index 1c23b9afb..99b2834f1 100644 --- a/include/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp @@ -209,10 +209,11 @@ struct traversal_ring_creator if (start_turn.is_clustered()) { - turn_type const& turn = m_turns[current_turn_index]; - if (turn.cluster_id == start_turn.cluster_id) + turn_type& turn = m_turns[current_turn_index]; + turn_operation_type& op = turn.operations[current_op_index]; + if (turn.cluster_id == start_turn.cluster_id + && op.enriched.get_next_turn_index() == start_turn_index) { - turn_operation_type& op = m_turns[start_turn_index].operations[current_op_index]; op.visited.set_finished(); m_visitor.visit_traverse(m_turns, m_turns[current_turn_index], start_op, "Early finish (cluster)"); return traverse_error_none; @@ -307,6 +308,23 @@ struct traversal_ring_creator } } + int get_operation_index(turn_type const& turn) const + { + // When starting with a continue operation, the one + // with the smallest (for intersection) or largest (for union) + // remaining distance (#8310b) + // Also to avoid skipping a turn in between, which can happen + // in rare cases (e.g. #130) + static const bool is_union + = operation_from_overlay::value == operation_union; + + turn_operation_type const& op0 = turn.operations[0]; + turn_operation_type const& op1 = turn.operations[1]; + return op0.remaining_distance <= op1.remaining_distance + ? (is_union ? 1 : 0) + : (is_union ? 0 : 1); + } + template void iterate(Rings& rings, std::size_t& finalized_ring_size, typename Backtrack::state_type& state) @@ -323,15 +341,8 @@ struct traversal_ring_creator if (turn.both(operation_continue)) { - // Traverse only one turn, the one with the SMALLEST remaining distance - // to avoid skipping a turn in between, which can happen in rare cases - // (e.g. #130) - turn_operation_type const& op0 = turn.operations[0]; - turn_operation_type const& op1 = turn.operations[1]; - int const op_index - = op0.remaining_distance <= op1.remaining_distance ? 0 : 1; - - traverse_with_operation(turn, turn_index, op_index, + traverse_with_operation(turn, turn_index, + get_operation_index(turn), rings, finalized_ring_size, state); } else @@ -374,13 +385,8 @@ struct traversal_ring_creator if (turn.both(operation_continue)) { - // Traverse only one turn, the one with the SMALLEST remaining distance - // to avoid skipping a turn in between, which can happen in rare cases - // (e.g. #130) - int const op_index - = op0.remaining_distance <= op1.remaining_distance ? 0 : 1; - - traverse_with_operation(turn, turn_index, op_index, + traverse_with_operation(turn, turn_index, + get_operation_index(turn), rings, finalized_ring_size, state); } else diff --git a/test/algorithms/overlay/multi_overlay_cases.hpp b/test/algorithms/overlay/multi_overlay_cases.hpp index 219957456..68215f842 100644 --- a/test/algorithms/overlay/multi_overlay_cases.hpp +++ b/test/algorithms/overlay/multi_overlay_cases.hpp @@ -1323,6 +1323,18 @@ static std::string pie_7_2_1_0_15[2] = "MULTIPOLYGON(((2500 2500,2791 3586,3062 3474,2500 2500)),((2500 2500,3474 3062,3586 2791,3625 2500,3586 2208,3474 1937,3295 1704,3062 1525,2791 1413,2499 1375,2208 1413,1937 1525,1704 1704,1525 1937,1413 2208,1375 2500,2500 2500)))" }; +static std::string case_precision_m1[2] = +{ + "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))", + "MULTIPOLYGON(((-1 -1,-1 8,2 8,2 7,2 3,4.0000005 2.9999995,4 7,4 8,8 8,8 -1,-1 -1)))" +}; + +static std::string case_precision_m2[2] = +{ + "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)),((3 6,3 7.5,4.5 7.5,4.5 6,3 6)))", + "MULTIPOLYGON(((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4.0000005 2.9999995,4 7,2 7)))" +}; + // Case, not literally on this list but derived, to mix polygon/multipolygon in call to difference static std::string ggl_list_20111025_vd[4] = { diff --git a/test/algorithms/overlay/overlay_cases.hpp b/test/algorithms/overlay/overlay_cases.hpp index c917af845..a84812a2d 100644 --- a/test/algorithms/overlay/overlay_cases.hpp +++ b/test/algorithms/overlay/overlay_cases.hpp @@ -700,6 +700,88 @@ static std::string dz_4[2] = { "POLYGON((20.486696243286133 60.650150299072266,24.282432556152344 49.304500579833984,34.362251281738281 55.748767852783203,30.764263153076172 44.3388671875,42.706855773925781 43.627620697021484,33.089447021484375 36.511661529541016,42.333145141601563 28.916570663452148,30.369846343994141 28.81260871887207,33.383872985839844 17.234743118286133,23.644252777099609 24.182485580444336,19.277351379394531 13.044195175170898,15.48161506652832 24.389842987060547,5.40179443359375 17.945577621459961,8.9997835159301758 29.355476379394531,-2.9428071975708008 30.06672477722168,6.6745977401733398 37.182682037353516,-2.5690991878509521 44.777774810791016,9.394200325012207 44.881736755371094,6.3801741600036621 56.459602355957031,16.119794845581055 49.511859893798828,20.486696243286133 60.650150299072266))" }; +static std::string case_precision_1[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4.000005 2.99999,2 3,2 7))" +}; + +static std::string case_precision_2[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4 2.999995,2 3,2 7))" +}; + +static std::string case_precision_3[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4.0000001 2.99999995,2 3,2 7))" +}; + +static std::string case_precision_4[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4 3.00000001,2 3,2 7))" +}; + +static std::string case_precision_5[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4 3,2.0000005 2.9999995,2 7))" +}; + +static std::string case_precision_6[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((-1 -1,-1 8,2 8,2 7,2 3,4.0000005 2.9999995,4 7,4 8,8 8,8 -1,-1 -1))" +}; + +static std::string case_precision_7[2] = +{ + // Needs larger margin for sectionalize (long double only) + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((2 7,4 7,4 3.00000002,2 3,2 7))" +}; + +static std::string case_precision_8[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4.00000006 3.00000009,4 7,2 7))" +}; + +static std::string case_precision_9[2] = +{ + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,3.99999 2.999995,4 7,2 7))" +}; + +static std::string case_precision_10[2] = +{ + // Needs 1.0e-5 for threshold in double + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4.000006 2.999991,4 7,2 7))" +}; + +static std::string case_precision_11[2] = +{ + // Needs ~0.5 for threshold in side_by_generic_form + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4.00000000000000089 2.99999999999999201,4 7,2 7))" +}; + +static std::string case_precision_12[2] = +{ + // Needs threshold for threshold a2 + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((1 1,2.99999999999999731e-12 1.00000000001,2.99999999999999731e-12 3.00000000001,1 3,1 1))" +}; + +static std::string case_precision_13[2] = +{ + // Needs threshold for threshold a2 + "POLYGON((0 0,0 4,2 4,2 3,4 3,4 0,0 0))", + "POLYGON((1 1,9.99999999999999912e-06 1,9.99999999999999912e-06 3,1 3,1 1))" +}; // ticket_17 is keyholed, so has a hole formed by an deliberate intersection // This will fail the intersection/traversal process diff --git a/test/algorithms/set_operations/difference/difference.cpp b/test/algorithms/set_operations/difference/difference.cpp index 51ff33a2b..289530e5f 100644 --- a/test/algorithms/set_operations/difference/difference.cpp +++ b/test/algorithms/set_operations/difference/difference.cpp @@ -304,15 +304,17 @@ void test_all() // Isovist - the # output polygons differ per compiler/pointtype, (very) small // rings might be discarded. We check area only + + // SQL Server gives: 0.279121891701124 and 224.889211358929 + // PostGIS gives: 0.279121991127244 and 224.889205853156 + // No robustness gives: 0.279121991127106 and 224.825363749290 + test_one("isovist", isovist1[0], isovist1[1], -1, -1, 0.279132, -1, -1, 224.8892, settings); } - // SQL Server gives: 0.279121891701124 and 224.889211358929 - // PostGIS gives: 0.279121991127244 and 224.889205853156 - // No robustness gives: 0.279121991127106 and 224.825363749290 #ifdef BOOST_GEOMETRY_TEST_INCLUDE_FAILING_TESTS test_one("geos_1", diff --git a/test/algorithms/set_operations/difference/difference_multi.cpp b/test/algorithms/set_operations/difference/difference_multi.cpp index 25e82b1c3..389216892 100644 --- a/test/algorithms/set_operations/difference/difference_multi.cpp +++ b/test/algorithms/set_operations/difference/difference_multi.cpp @@ -463,6 +463,11 @@ void test_areal() TEST_DIFFERENCE(case_recursive_boxes_87, 4, 2.0, 4, 2.5, 8); TEST_DIFFERENCE(case_recursive_boxes_88, 3, 4.75, 5, 6.75, 4); + // Output of A can be 0 or 1 polygons (with a very small area) + TEST_DIFFERENCE(case_precision_m1, -1, 0.0, 1, 57.0, -1); + // Output of A can be 1 or 2 polygons (one with a very small area) + TEST_DIFFERENCE(case_precision_m2, -1, 1.0, 1, 57.75, -1); + { ut_settings sym_settings; #if defined(BOOST_GEOMETRY_NO_ROBUSTNESS) diff --git a/test/algorithms/set_operations/difference/test_difference.hpp b/test/algorithms/set_operations/difference/test_difference.hpp index 3dc31c7e9..4a505bcdb 100644 --- a/test/algorithms/set_operations/difference/test_difference.hpp +++ b/test/algorithms/set_operations/difference/test_difference.hpp @@ -273,7 +273,15 @@ std::string test_difference(std::string const& caseid, G1 const& g1, G2 const& g ); } - BOOST_CHECK_CLOSE(area, expected_area, settings.percentage); + if (expected_area > 0) + { + BOOST_CHECK_CLOSE(area, expected_area, settings.percentage); + } + else + { + // Compare 0 with 0 or a very small detected area + BOOST_CHECK_LE(area, settings.percentage); + } #endif diff --git a/test/algorithms/set_operations/intersection/intersection.cpp b/test/algorithms/set_operations/intersection/intersection.cpp index f0934cb98..c59c4afff 100644 --- a/test/algorithms/set_operations/intersection/intersection.cpp +++ b/test/algorithms/set_operations/intersection/intersection.cpp @@ -46,6 +46,16 @@ BOOST_GEOMETRY_REGISTER_LINESTRING_TEMPLATED(std::vector) (test_one) \ ( #caseid, caseid[0], caseid[1], clips, points, area) +#define TEST_INTERSECTION_REV(caseid, clips, points, area) \ + (test_one) \ + ( #caseid "_rev", caseid[1], caseid[0], clips, points, area) + +#define TEST_INTERSECTION_WITH(caseid, index1, index2, \ + clips, points, area, settings) \ + (test_one) \ + ( #caseid #index1 "_" #index2, caseid[index1], caseid[index2], \ + clips, points, area, settings) + #if defined(BOOST_GEOMETRY_NO_SELF_TURNS) #define TEST_INTERSECTION_IGNORE(caseid, clips, points, area) \ { ut_settings ignore_validity; ignore_validity.test_validity = false; \ @@ -197,8 +207,7 @@ void test_areal() // In most cases: 0 (no intersection) // In some cases: 1.430511474609375e-05 (clang/gcc on Xubuntu using b2) // In some cases: 5.6022983000000002e-05 (powerpc64le-gcc-6-0) - test_one("geos_2", - geos_2[0], geos_2[1], + test_one("geos_2", geos_2[0], geos_2[1], 0, 0, 6.0e-5, ut_settings(-1.0)); // -1 denotes: compare with <= #if ! defined(BOOST_GEOMETRY_NO_ROBUSTNESS) @@ -276,12 +285,12 @@ void test_areal() test_one("ticket_8652", ticket_8652[0], ticket_8652[1], 1, 4, 0.0003); - test_one("ticket_8310a", ticket_8310a[0], ticket_8310a[1], - 1, 5, 0.3843747); - test_one("ticket_8310b", ticket_8310b[0], ticket_8310b[1], - 1, 5, 0.3734379); - test_one("ticket_8310c", ticket_8310c[0], ticket_8310c[1], - 1, 5, 0.4689541); + TEST_INTERSECTION(ticket_8310a, 1, 5, 0.3843747); + TEST_INTERSECTION(ticket_8310b, 1, 5, 0.3734379); + TEST_INTERSECTION(ticket_8310c, 1, 5, 0.4689541); + TEST_INTERSECTION_REV(ticket_8310a, 1, 5, 0.3843747); + TEST_INTERSECTION_REV(ticket_8310b, 1, 5, 0.3734379); + TEST_INTERSECTION_REV(ticket_8310c, 1, 5, 0.4689541); test_one("ticket_9081_15", ticket_9081_15[0], ticket_9081_15[1], @@ -296,7 +305,7 @@ void test_areal() // mingw 5.6022954e-5 test_one("ticket_10108_b", ticket_10108_b[0], ticket_10108_b[1], - 0, 0, 5.6022983e-5); + 0, 0, 5.6022983e-5, ut_settings(-1.0)); #endif test_one("ticket_10747_a", @@ -355,6 +364,37 @@ void test_areal() TEST_INTERSECTION(case_105, 1, 34, 76.0); + TEST_INTERSECTION(case_precision_1, 0, 0, 0.0); + TEST_INTERSECTION(case_precision_2, 0, 0, 0.0); + TEST_INTERSECTION(case_precision_3, 0, 0, 0.0); + TEST_INTERSECTION(case_precision_4, 0, 0, 0.0); + TEST_INTERSECTION(case_precision_5, 0, 0, 0.0); + TEST_INTERSECTION(case_precision_6, 1, -1, 14.0); + TEST_INTERSECTION(case_precision_7, 0, -1, 0.0); + TEST_INTERSECTION(case_precision_8, 1, -1, 14.0); + TEST_INTERSECTION(case_precision_9, 1, -1, 14.0); + TEST_INTERSECTION(case_precision_10, 1, -1, 14.0); + TEST_INTERSECTION(case_precision_11, 1, -1, 14.0); + + TEST_INTERSECTION_REV(case_precision_1, 0, 0, 0.0); + TEST_INTERSECTION_REV(case_precision_2, 0, 0, 0.0); + TEST_INTERSECTION_REV(case_precision_3, 0, 0, 0.0); + TEST_INTERSECTION_REV(case_precision_4, 0, 0, 0.0); + TEST_INTERSECTION_REV(case_precision_5, 0, 0, 0.0); + TEST_INTERSECTION_REV(case_precision_6, 1, -1, 14.0); + TEST_INTERSECTION_REV(case_precision_7, 0, -1, 0.0); + TEST_INTERSECTION_REV(case_precision_8, 1, -1, 14.0); + TEST_INTERSECTION_REV(case_precision_9, 1, -1, 14.0); + TEST_INTERSECTION_REV(case_precision_10, 1, -1, 14.0); + TEST_INTERSECTION_REV(case_precision_11, 1, -1, 14.0); + { + ut_settings settings(0.01); + TEST_INTERSECTION_WITH(case_precision_12, 0, 1, 1, -1, 2.0, settings); + TEST_INTERSECTION_WITH(case_precision_13, 0, 1, 1, -1, 2.0, settings); + TEST_INTERSECTION_WITH(case_precision_12, 1, 0, 1, -1, 2.0, settings); + TEST_INTERSECTION_WITH(case_precision_13, 1, 0, 1, -1, 2.0, settings); + } + #ifndef BOOST_GEOMETRY_NO_SELF_TURNS TEST_INTERSECTION(case_106, 2, -1, 3.5); TEST_INTERSECTION(case_107, 3, -1, 3.0); @@ -890,7 +930,6 @@ int test_main(int, char* []) test_all >(); #endif -#endif // Commented, because exception is now disabled: // test_exception >(); @@ -918,6 +957,7 @@ int test_main(int, char* []) #if defined(BOOST_HAS_LONG_LONG) test_ticket_10868("MULTIPOLYGON(((33520458 6878575,33480192 14931538,31446819 18947953,30772384 19615678,30101303 19612322,30114725 16928001,33520458 6878575)))"); #endif +#endif #endif return 0; diff --git a/test/algorithms/set_operations/intersection/intersection_multi.cpp b/test/algorithms/set_operations/intersection/intersection_multi.cpp index a0f4a215f..139507780 100644 --- a/test/algorithms/set_operations/intersection/intersection_multi.cpp +++ b/test/algorithms/set_operations/intersection/intersection_multi.cpp @@ -32,6 +32,10 @@ (test_one) \ ( #caseid, caseid[0], caseid[1], clips, points, area) +#define TEST_INTERSECTION_REV(caseid, clips, points, area) \ + (test_one) \ + ( #caseid "_rev", caseid[1], caseid[0], clips, points, area) + #define TEST_INTERSECTION_IGNORE(caseid, clips, points, area) \ { ut_settings ignore_validity; ignore_validity.test_validity = false; \ (test_one) \ @@ -368,6 +372,17 @@ void test_areal() TEST_INTERSECTION(case_recursive_boxes_87, 0, -1, 0.0); TEST_INTERSECTION(case_recursive_boxes_88, 4, -1, 3.5); +#ifndef BOOST_GEOMETRY_NO_ROBUSTNESS + TEST_INTERSECTION(case_precision_m1, 1, -1, 14.0); + TEST_INTERSECTION(case_precision_m2, 2, -1, 15.25); + TEST_INTERSECTION_REV(case_precision_m1, 1, -1, 14.0); + TEST_INTERSECTION_REV(case_precision_m2, 2, -1, 15.25); +#else + // Validity: false positives (very small triangles looking like a line) + TEST_INTERSECTION_IGNORE(case_precision_m1, 1, -1, 14.0); + TEST_INTERSECTION_IGNORE(case_precision_m2, 2, -1, 15.25); +#endif + test_one("ggl_list_20120915_h2_a", ggl_list_20120915_h2[0], ggl_list_20120915_h2[1], 2, 10, 6.0); // Area from SQL Server diff --git a/test/algorithms/set_operations/intersection/test_intersection.hpp b/test/algorithms/set_operations/intersection/test_intersection.hpp index 8a9d6ac88..d780ee14e 100644 --- a/test/algorithms/set_operations/intersection/test_intersection.hpp +++ b/test/algorithms/set_operations/intersection/test_intersection.hpp @@ -151,7 +151,15 @@ check_result( double const detected_length_or_area = boost::numeric_cast(length_or_area); if (settings.percentage > 0.0) { - BOOST_CHECK_CLOSE(detected_length_or_area, expected_length_or_area, settings.percentage); + if (expected_length_or_area > 0) + { + BOOST_CHECK_CLOSE(detected_length_or_area, expected_length_or_area, settings.percentage); + } + else + { + // Compare 0 with 0 or a very small detected area + BOOST_CHECK_LE(detected_length_or_area, settings.percentage); + } } else { diff --git a/test/algorithms/set_operations/union/union.cpp b/test/algorithms/set_operations/union/union.cpp index e989e6187..d2dee0662 100644 --- a/test/algorithms/set_operations/union/union.cpp +++ b/test/algorithms/set_operations/union/union.cpp @@ -26,6 +26,14 @@ #include +#define TEST_UNION(caseid, clips, holes, points, area) \ + (test_one) \ + ( #caseid, caseid[0], caseid[1], clips, holes, points, area) + +#define TEST_UNION_REV(caseid, clips, holes, points, area) \ + (test_one) \ + ( #caseid "_rev", caseid[1], caseid[0], clips, holes, points, area) + template void test_areal() @@ -247,6 +255,34 @@ void test_areal() test_one("108", case_108[0], case_108[1], 1, 0, 13, 5.0); + TEST_UNION(case_precision_1, 1, 0, -1, 22.0); + TEST_UNION(case_precision_2, 1, 0, -1, 22.0); + TEST_UNION(case_precision_3, 1, 0, -1, 22.0); + TEST_UNION(case_precision_4, 1, 0, -1, 22.0); + TEST_UNION(case_precision_5, 1, 0, -1, 22.0); + TEST_UNION(case_precision_6, 1, 0, -1, 71.0); + TEST_UNION(case_precision_7, 1, 0, -1, 22.0); + TEST_UNION(case_precision_8, 1, 1, -1, 73.0); + TEST_UNION(case_precision_9, 1, 1, -1, 73.0); + TEST_UNION(case_precision_10, 1, 1, -1, 73.0); + TEST_UNION(case_precision_11, 1, 1, -1, 73.0); + TEST_UNION(case_precision_12, 1, 0, -1, 14.0); + TEST_UNION(case_precision_13, 1, 0, -1, 14.0); + + TEST_UNION_REV(case_precision_1, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_2, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_3, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_4, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_5, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_6, 1, 0, -1, 71.0); + TEST_UNION_REV(case_precision_7, 1, 0, -1, 22.0); + TEST_UNION_REV(case_precision_8, 1, 1, -1, 73.0); + TEST_UNION_REV(case_precision_9, 1, 1, -1, 73.0); + TEST_UNION_REV(case_precision_10, 1, 1, -1, 73.0); + TEST_UNION_REV(case_precision_11, 1, 1, -1, 73.0); + TEST_UNION_REV(case_precision_12, 1, 0, -1, 14.0); + TEST_UNION_REV(case_precision_13, 1, 0, -1, 14.0); + /* test_one(102, simplex_normal[0], simplex_reversed[1], @@ -342,22 +378,23 @@ void test_areal() ticket_5103[0], ticket_5103[1], 1, 0, 25, 2515271327070.5); - test_one("ticket_8310a", ticket_8310a[0], ticket_8310a[1], - 1, 0, 5, 10.5000019595); - test_one("ticket_8310b", ticket_8310b[0], ticket_8310b[1], - 1, 0, 5, 10.5000019595); - test_one("ticket_8310c", ticket_8310c[0], ticket_8310c[1], - 1, 0, 5, 10.5000019595); + TEST_UNION(ticket_8310a, 1, 0, 5, 10.5000019595); + TEST_UNION(ticket_8310b, 1, 0, 5, 10.5000019595); + TEST_UNION(ticket_8310c, 1, 0, 5, 10.5000019595); + TEST_UNION_REV(ticket_8310a, 1, 0, 5, 10.5000019595); + TEST_UNION_REV(ticket_8310b, 1, 0, 5, 10.5000019595); + TEST_UNION_REV(ticket_8310c, 1, 0, 5, 10.5000019595); test_one("ticket_9081_15", ticket_9081_15[0], ticket_9081_15[1], - 1, 0, 10, 0.0403425433); + 1, 0, -1, 0.0403425433); test_one("ticket_9563", ticket_9563[0], ticket_9563[1], 1, 0, 13, 150.0); + // Float result is OK but a bit larger test_one("ticket_9756", ticket_9756[0], ticket_9756[1], - 1, 0, 10, 1289.08374); + 1, 0, 10, if_typed(1291.5469, 1289.08374)); #if defined(BOOST_GEOMETRY_NO_ROBUSTNESS) test_one("ticket_10108_a", ticket_10108_a[0], ticket_10108_a[1], @@ -373,7 +410,7 @@ void test_areal() #endif test_one("ticket_10866", ticket_10866[0], ticket_10866[1], - 1, 0, 14, 332760303.5); + 1, 0, 14, if_typed(332752493.0, 332760303.5)); test_one("ticket_11725", ticket_11725[0], ticket_11725[1], 1, 1, 10, 7.5); @@ -398,17 +435,17 @@ void test_areal() // Robustness issues, followed out buffer-robustness-tests, test them also reverse #if ! defined(BOOST_GEOMETRY_NO_ROBUSTNESS) test_one("buffer_rt_f", buffer_rt_f[0], buffer_rt_f[1], - 1, 0, 15, 4.60853); + 1, 0, -1, 4.60853); test_one("buffer_rt_f_rev", buffer_rt_f[1], buffer_rt_f[0], - 1, 0, 15, 4.60853); + 1, 0, -1, 4.60853); test_one("buffer_rt_g", buffer_rt_g[0], buffer_rt_g[1], - 1, 0, if_typed(16, 11), 16.571); + 1, 0, -1, 16.571); test_one("buffer_rt_g_rev", buffer_rt_g[1], buffer_rt_g[0], - 1, 0, if_typed(16, 11), 16.571); + 1, 0, -1, 16.571); test_one("buffer_rt_i", buffer_rt_i[0], buffer_rt_i[1], - 1, 0, if_typed(11, 13), 13.6569); + 1, 0, -1, 13.6569); test_one("buffer_rt_i_rev", buffer_rt_i[1], buffer_rt_i[0], - 1, 0, 13, 13.6569); + 1, 0, -1, 13.6569); #endif test_one("buffer_rt_j", buffer_rt_j[0], buffer_rt_j[1], @@ -429,23 +466,23 @@ void test_areal() 1, 0, 9, 19.4852); test_one("buffer_rt_m2", buffer_rt_m2[0], buffer_rt_m2[1], - 1, 0, 12, 21.4853); + 1, 0, -1, 21.4853); test_one("buffer_rt_m2_rev", buffer_rt_m2[1], buffer_rt_m2[0], 1, 0, 15, 21.4853); #if ! defined(BOOST_GEOMETRY_NO_ROBUSTNESS) test_one("buffer_rt_q", buffer_rt_q[0], buffer_rt_q[1], - 1, 0, if_typed(16, 12), 18.5710); + 1, 0, -1, 18.5710); test_one("buffer_rt_q_rev", buffer_rt_q[1], buffer_rt_q[0], - 1, 0, if_typed(16, 12), 18.5710); + 1, 0, -1, 18.5710); test_one("buffer_rt_r", buffer_rt_r[0], buffer_rt_r[1], - 1, 0, if_typed(18, 14), 21.07612); + 1, 0, -1, 21.07612); test_one("buffer_rt_r_rev", buffer_rt_r[1], buffer_rt_r[0], - 1, 0, if_typed(18, 14), 21.07612); + 1, 0, -1, 21.07612); test_one("buffer_rt_t", buffer_rt_t[0], buffer_rt_t[1], - 1, 0, 9, 15.6569); + 1, 0, -1, 15.6569); test_one("buffer_rt_t_rev", buffer_rt_t[1], buffer_rt_t[0], - 1, 0, 10, 15.6569); + 1, 0, -1, 15.6569); #endif test_one("buffer_mp1", buffer_mp1[0], buffer_mp1[1],