From ec7f9c9887ad4b8a41f09b7d392c98aafdc34321 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Wed, 15 Feb 2023 10:48:03 +0100 Subject: [PATCH] [union] fix missing interior ring and double traversed exterior ring fixes: #1109 and #1108 keeps fixed: #1081 --- .../detail/overlay/colocate_clusters.hpp | 107 ++++++++++++++++++ .../detail/overlay/get_clusters.hpp | 5 +- .../detail/overlay/handle_colocations.hpp | 26 ++--- .../algorithms/detail/overlay/traversal.hpp | 48 +++++--- .../buffer/buffer_multi_polygon.cpp | 11 +- test/algorithms/overlay/get_clusters.cpp | 12 +- .../overlay/multi_overlay_cases.hpp | 6 + test/algorithms/overlay/overlay_cases.hpp | 6 + .../algorithms/set_operations/union/union.cpp | 6 + .../set_operations/union/union_multi.cpp | 2 + 10 files changed, 183 insertions(+), 46 deletions(-) create mode 100644 include/boost/geometry/algorithms/detail/overlay/colocate_clusters.hpp diff --git a/include/boost/geometry/algorithms/detail/overlay/colocate_clusters.hpp b/include/boost/geometry/algorithms/detail/overlay/colocate_clusters.hpp new file mode 100644 index 000000000..06d6e7bd6 --- /dev/null +++ b/include/boost/geometry/algorithms/detail/overlay/colocate_clusters.hpp @@ -0,0 +1,107 @@ +// Boost.Geometry + +// Copyright (c) 2023 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP + +#include +#include +#include +#include + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace overlay +{ + +// Default implementation, using the first point for all turns in the cluster. +template +< + typename Point, + typename CoordinateType = typename geometry::coordinate_type::type, + typename CsTag = typename geometry::cs_tag::type, + bool IsIntegral = std::is_integral::value +> +struct cluster_colocator +{ + template + static inline void apply(TurnIndices const& indices, Turns& turns) + { + // This approach works for all but one testcase (rt_p13) + // The problem is fill_sbs, which uses sides and these sides might change slightly + // depending on the exact location of the cluster. + // Using the centroid is, on the average, a safer choice for sides. + // Alternatively fill_sbs could be revised, but that requires a lot of work + // and is outside current scope. + // Integer coordinates are always colocated already and do not need centroid calculation. + // Geographic/spherical coordinates might (in extremely rare cases) cross the date line + // and therefore the first point is taken for them as well. + auto it = indices.begin(); + auto const& first_point = turns[*it].point; + for (++it; it != indices.end(); ++it) + { + turns[*it].point = first_point; + } + } +}; + +// Specialization for non-integral cartesian coordinates, calculating +// the centroid of the points of the turns in the cluster. +template +struct cluster_colocator +{ + template + static inline void apply(TurnIndices const& indices, Turns& turns) + { + CoordinateType centroid_0 = 0; + CoordinateType centroid_1 = 0; + for (const auto& index : indices) + { + centroid_0 += geometry::get<0>(turns[index].point); + centroid_1 += geometry::get<1>(turns[index].point); + } + centroid_0 /= indices.size(); + centroid_1 /= indices.size(); + for (const auto& index : indices) + { + geometry::set<0>(turns[index].point, centroid_0); + geometry::set<1>(turns[index].point, centroid_1); + } + } +}; + +// Moves intersection points per cluster such that they are identical. +// Because clusters are intersection close together, and +// handled as one location. Then they should also have one location. +// It is necessary to avoid artefacts and invalidities. +template +inline void colocate_clusters(Clusters const& clusters, Turns& turns) +{ + for (auto const& pair : clusters) + { + auto const& indices = pair.second.turn_indices; + if (indices.size() < 2) + { + // Defensive check + continue; + } + using point_t = decltype(turns[*indices.begin()].point); + cluster_colocator::apply(indices, turns); + } +} + + +}} // namespace detail::overlay +#endif //DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP diff --git a/include/boost/geometry/algorithms/detail/overlay/get_clusters.hpp b/include/boost/geometry/algorithms/detail/overlay/get_clusters.hpp index 9e2627f36..2747fa68b 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_clusters.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_clusters.hpp @@ -39,9 +39,8 @@ private: template static inline T threshold() { - // Tuned by cases of #1081. It should just be one epsilon to distinguish between - // different turn-points and real colocated clusters. - return T(1); + // Points within some epsilons are considered as equal. + return T(100); } public: // Returns true if point are considered equal (within an epsilon) diff --git a/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp index 7803db970..ac64c22cf 100644 --- a/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -52,22 +53,22 @@ namespace boost { namespace geometry namespace detail { namespace overlay { +// Removes clusters which have only one point left, or are empty. template inline void remove_clusters(Turns& turns, Clusters& clusters) { - typename Clusters::iterator it = clusters.begin(); + auto it = clusters.begin(); while (it != clusters.end()) { // Hold iterator and increase. We can erase cit, this keeps the // iterator valid (cf The standard associative-container erase idiom) - typename Clusters::iterator current_it = it; + auto current_it = it; ++it; - std::set const& turn_indices - = current_it->second.turn_indices; + auto const& turn_indices = current_it->second.turn_indices; if (turn_indices.size() == 1) { - signed_size_type const turn_index = *turn_indices.begin(); + auto const turn_index = *turn_indices.begin(); turns[turn_index].cluster_id = -1; clusters.erase(current_it); } @@ -78,18 +79,16 @@ template inline void cleanup_clusters(Turns& turns, Clusters& clusters) { // Removes discarded turns from clusters - for (typename Clusters::iterator mit = clusters.begin(); - mit != clusters.end(); ++mit) + for (auto& pair : clusters) { - cluster_info& cinfo = mit->second; - std::set& ids = cinfo.turn_indices; - for (std::set::iterator sit = ids.begin(); - sit != ids.end(); /* no increment */) + auto& cinfo = pair.second; + auto& ids = cinfo.turn_indices; + for (auto sit = ids.begin(); sit != ids.end(); /* no increment */) { - std::set::iterator current_it = sit; + auto current_it = sit; ++sit; - signed_size_type const turn_index = *current_it; + auto const turn_index = *current_it; if (turns[turn_index].discarded) { ids.erase(current_it); @@ -98,6 +97,7 @@ inline void cleanup_clusters(Turns& turns, Clusters& clusters) } remove_clusters(turns, clusters); + colocate_clusters(clusters, turns); } template diff --git a/include/boost/geometry/algorithms/detail/overlay/traversal.hpp b/include/boost/geometry/algorithms/detail/overlay/traversal.hpp index b65f00424..1e1dd6739 100644 --- a/include/boost/geometry/algorithms/detail/overlay/traversal.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/traversal.hpp @@ -968,31 +968,43 @@ public : { turn_type const& current_turn = m_turns[turn_index]; + bool const back_at_start_cluster + = has_points + && current_turn.is_clustered() + && m_turns[start_turn_index].cluster_id == current_turn.cluster_id; if (BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection)) { - if (has_points) - { - bool const back_at_start_cluster - = current_turn.is_clustered() - && m_turns[start_turn_index].cluster_id == current_turn.cluster_id; + // Intersection or difference - if (turn_index == start_turn_index || back_at_start_cluster) - { - // Intersection can always be finished if returning - turn_index = start_turn_index; - op_index = start_op_index; - return true; - } + if (has_points && (turn_index == start_turn_index || back_at_start_cluster)) + { + // Intersection can always be finished if returning + turn_index = start_turn_index; + op_index = start_op_index; + return true; } if (! current_turn.is_clustered() - && current_turn.both(operation_intersection)) - { - if (analyze_ii_intersection(turn_index, op_index, + && current_turn.both(operation_intersection) + && analyze_ii_intersection(turn_index, op_index, current_turn, previous_seg_id)) - { - return true; - } + { + return true; + } + } + else if (turn_index == start_turn_index || back_at_start_cluster) + { + // Union or buffer: cannot return immediately to starting turn, because it then + // might miss a formed multi polygon with a touching point. + auto const& current_op = current_turn.operations[op_index]; + signed_size_type const next_turn_index = current_op.enriched.get_next_turn_index(); + bool const to_other_turn = next_turn_index >= 0 && m_turns[next_turn_index].cluster_id != current_turn.cluster_id; + if (! to_other_turn) + { + // Return to starting point + turn_index = start_turn_index; + op_index = start_op_index; + return true; } } diff --git a/test/algorithms/buffer/buffer_multi_polygon.cpp b/test/algorithms/buffer/buffer_multi_polygon.cpp index 1ce3c9e8b..64f8691a0 100644 --- a/test/algorithms/buffer/buffer_multi_polygon.cpp +++ b/test/algorithms/buffer/buffer_multi_polygon.cpp @@ -540,7 +540,10 @@ void test_all() test_one("rt_p11", rt_p11, join_miter, end_flat, 28.7426, 1.0); test_one("rt_p12", rt_p12, join_miter, end_flat, 22.5711, 1.0); + + // Needs centroid of cluster turn points test_one("rt_p13", rt_p13, join_miter, end_flat, 19.9142, 1.0); + test_one("rt_p14", rt_p14, join_miter, end_flat, 20.8284, 1.0); test_one("rt_p15", rt_p15, join_miter, end_flat, 23.6569, 1.0); test_one("rt_p16", rt_p16, join_miter, end_flat, 23.4853, 1.0); @@ -631,11 +634,9 @@ void test_all() 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) - // Fails since get_cluster works with an epsilon of 1 (instead of 1000 before). - // With 3 it still succeeds but that causes regression in issue #issue_1081b + // Needs an espilon in get_cluster of 3 or higher test_one("nores_1ea1", nores_1ea1, join_round32, end_flat, 28.9755, 1.0); -#endif + test_one("nores_804e", nores_804e, join_round32, end_flat, 26.4503, 1.0); test_one("nores_51c6", nores_51c6, join_round32, end_flat, 20.2419, 1.0); test_one("nores_e5f3", nores_e5f3, join_round32, end_flat, 14.5503, 1.0); @@ -651,7 +652,7 @@ void test_all() #if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES) // Erroneous case with rescaling test_one("res_8618", res_8618, join_round32, end_flat, 48.1085, 1.0); - #endif +#endif test_one("res_3b4d", res_3b4d, join_round32, end_flat, 48.4739, 1.0); test_one("neighbouring_small", diff --git a/test/algorithms/overlay/get_clusters.cpp b/test/algorithms/overlay/get_clusters.cpp index 7cb275f11..5b9a3602f 100644 --- a/test/algorithms/overlay/get_clusters.cpp +++ b/test/algorithms/overlay/get_clusters.cpp @@ -101,13 +101,11 @@ int test_main(int, char* []) test_get_clusters(); test_get_clusters(); - // This constant relates to the threshold in get_clusters, - // which is now return T(1) (earlier it was 1000) - double const multiplier = 1.0 / 1000.0; - - test_get_clusters_border_cases(1.0e-4 * multiplier); - test_get_clusters_border_cases(1.0e-13 * multiplier); - test_get_clusters_border_cases(1.0e-16 * multiplier); + // These constant relate to the threshold in get_clusters.hpp, + // and the used floating point type. + test_get_clusters_border_cases(1.0e-5); + test_get_clusters_border_cases(1.0e-14); + test_get_clusters_border_cases(1.0e-17); return 0; } diff --git a/test/algorithms/overlay/multi_overlay_cases.hpp b/test/algorithms/overlay/multi_overlay_cases.hpp index 45fc7b532..d0a53ab9c 100644 --- a/test/algorithms/overlay/multi_overlay_cases.hpp +++ b/test/algorithms/overlay/multi_overlay_cases.hpp @@ -1574,6 +1574,12 @@ static std::string issue_930[2] = "MULTIPOLYGON(((-10 2,5 3,20 2,-10 2)))" }; +static std::string issue_1109[2] = +{ + "MULTIPOLYGON(((-63 -88,0 -115.40000152587890625,0 -124,-63 -124,-63 -88)),((0 -150,0 -124,13 -124,13 -121,0 -115.40000152587890625,0 -88,40 -88,40 -150,0 -150)))", + "MULTIPOLYGON(((0 -88,0 -115.40000152587890625,-10 -88,0 -88)))" +}; + static std::string bug_21155501[2] = { "MULTIPOLYGON(((-8.3935546875 27.449790329784214,4.9658203125 18.729501999072138,11.8212890625 23.563987128451217,9.7119140625 25.48295117535531,9.8876953125 31.728167146023935,8.3056640625 32.99023555965106,8.5693359375 37.16031654673677,-1.8896484375 35.60371874069731,-0.5712890625 32.02670629333614,-8.9208984375 29.458731185355344,-8.3935546875 27.449790329784214)))", diff --git a/test/algorithms/overlay/overlay_cases.hpp b/test/algorithms/overlay/overlay_cases.hpp index fc77f93e5..e9b9d1960 100644 --- a/test/algorithms/overlay/overlay_cases.hpp +++ b/test/algorithms/overlay/overlay_cases.hpp @@ -1101,6 +1101,12 @@ static std::string issue_1081c[2] = "POLYGON((423.217316892426425 67.5751428875900331,414.40388971336813 96.7777620478938587,410.429206867910466 100.02249750988706,412.848579034710781 101.010025694314706,412.502954439453561 102.98508206238165,416.996074177797027 109.05132662251431,427.364812035512671 72.6538592632950468,423.217316892426425 67.5751428875900331))" }; +static std::string issue_1108[2] = +{ + "POLYGON((17 -2,18 -1.269679093926235014,12 0,22 0,17 -2))", + "POLYGON((22 1,22 0,14 0,18 -1.2696790939262529996,12 0,22 1))" +}; + static std::string ggl_list_20120229_volker[3] = { "POLYGON((1716 1554,2076 2250,2436 2352,2796 1248,3156 2484,3516 2688,3516 2688,3156 2484,2796 1248,2436 2352,2076 2250, 1716 1554))", diff --git a/test/algorithms/set_operations/union/union.cpp b/test/algorithms/set_operations/union/union.cpp index 736b41bf0..623f59beb 100644 --- a/test/algorithms/set_operations/union/union.cpp +++ b/test/algorithms/set_operations/union/union.cpp @@ -416,7 +416,10 @@ void test_areal() TEST_UNION(ticket_10108_a, count_set(1, 2), 0, 8, 0.0435229); TEST_UNION(ticket_10108_b, count_set(1, 2), 0, 10, 2424.3449); +#if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES) + // With rescaling, there is a dependency on cluster tolerance, which alters the result. TEST_UNION(ticket_10866, 1, 0, 14, 332760303.5); +#endif TEST_UNION(ticket_11725, 1, 1, 10, 7.5); @@ -456,6 +459,9 @@ void test_areal() TEST_UNION(issue_1081c, 1, 1, -1, 2338.08); TEST_UNION_REV(issue_1081c, 1, 1, -1, 2338.08); + TEST_UNION(issue_1108, 1, 0, -1, 12.1742); + TEST_UNION_REV(issue_1108, 1, 0, -1, 12.1742); + { // Rescaling produces an invalid result ut_settings settings; diff --git a/test/algorithms/set_operations/union/union_multi.cpp b/test/algorithms/set_operations/union/union_multi.cpp index 29f421ebb..cc2f9309a 100644 --- a/test/algorithms/set_operations/union/union_multi.cpp +++ b/test/algorithms/set_operations/union/union_multi.cpp @@ -445,6 +445,8 @@ void test_areal() TEST_UNION(issue_888_34, 15, 0, -1, 0.3017459); TEST_UNION(issue_888_37, 52, 3, -1, 0.4033294); + TEST_UNION(issue_1109, 2, 0, -1, 3946.5); + // One or two polygons, the ideal case is 1 TEST_UNION(mail_2019_01_21_johan, count_set(1, 2), 0, -1, 0.00058896);