Skip to content

Commit dfcbb75

Browse files
authored
Merge pull request #1113 from barendgehrels/issue-1108-fix-interior-ring-for-union
[union] fix missing interior ring and double traversed exterior ring
2 parents ee3ef77 + ec7f9c9 commit dfcbb75

File tree

10 files changed

+183
-46
lines changed

10 files changed

+183
-46
lines changed
Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
// Boost.Geometry
2+
3+
// Copyright (c) 2023 Barend Gehrels, Amsterdam, the Netherlands.
4+
5+
// Use, modification and distribution is subject to the Boost Software License,
6+
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7+
// http://www.boost.org/LICENSE_1_0.txt)
8+
9+
#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP
10+
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP
11+
12+
#include <boost/geometry/core/access.hpp>
13+
#include <boost/geometry/core/cs.hpp>
14+
#include <boost/geometry/core/coordinate_type.hpp>
15+
#include <boost/geometry/core/tags.hpp>
16+
17+
namespace boost { namespace geometry
18+
{
19+
20+
#ifndef DOXYGEN_NO_DETAIL
21+
namespace detail { namespace overlay
22+
{
23+
24+
// Default implementation, using the first point for all turns in the cluster.
25+
template
26+
<
27+
typename Point,
28+
typename CoordinateType = typename geometry::coordinate_type<Point>::type,
29+
typename CsTag = typename geometry::cs_tag<Point>::type,
30+
bool IsIntegral = std::is_integral<CoordinateType>::value
31+
>
32+
struct cluster_colocator
33+
{
34+
template <typename TurnIndices, typename Turns>
35+
static inline void apply(TurnIndices const& indices, Turns& turns)
36+
{
37+
// This approach works for all but one testcase (rt_p13)
38+
// The problem is fill_sbs, which uses sides and these sides might change slightly
39+
// depending on the exact location of the cluster.
40+
// Using the centroid is, on the average, a safer choice for sides.
41+
// Alternatively fill_sbs could be revised, but that requires a lot of work
42+
// and is outside current scope.
43+
// Integer coordinates are always colocated already and do not need centroid calculation.
44+
// Geographic/spherical coordinates might (in extremely rare cases) cross the date line
45+
// and therefore the first point is taken for them as well.
46+
auto it = indices.begin();
47+
auto const& first_point = turns[*it].point;
48+
for (++it; it != indices.end(); ++it)
49+
{
50+
turns[*it].point = first_point;
51+
}
52+
}
53+
};
54+
55+
// Specialization for non-integral cartesian coordinates, calculating
56+
// the centroid of the points of the turns in the cluster.
57+
template <typename Point, typename CoordinateType>
58+
struct cluster_colocator<Point, CoordinateType, geometry::cartesian_tag, false>
59+
{
60+
template <typename TurnIndices, typename Turns>
61+
static inline void apply(TurnIndices const& indices, Turns& turns)
62+
{
63+
CoordinateType centroid_0 = 0;
64+
CoordinateType centroid_1 = 0;
65+
for (const auto& index : indices)
66+
{
67+
centroid_0 += geometry::get<0>(turns[index].point);
68+
centroid_1 += geometry::get<1>(turns[index].point);
69+
}
70+
centroid_0 /= indices.size();
71+
centroid_1 /= indices.size();
72+
for (const auto& index : indices)
73+
{
74+
geometry::set<0>(turns[index].point, centroid_0);
75+
geometry::set<1>(turns[index].point, centroid_1);
76+
}
77+
}
78+
};
79+
80+
// Moves intersection points per cluster such that they are identical.
81+
// Because clusters are intersection close together, and
82+
// handled as one location. Then they should also have one location.
83+
// It is necessary to avoid artefacts and invalidities.
84+
template <typename Clusters, typename Turns>
85+
inline void colocate_clusters(Clusters const& clusters, Turns& turns)
86+
{
87+
for (auto const& pair : clusters)
88+
{
89+
auto const& indices = pair.second.turn_indices;
90+
if (indices.size() < 2)
91+
{
92+
// Defensive check
93+
continue;
94+
}
95+
using point_t = decltype(turns[*indices.begin()].point);
96+
cluster_colocator<point_t>::apply(indices, turns);
97+
}
98+
}
99+
100+
101+
}} // namespace detail::overlay
102+
#endif //DOXYGEN_NO_DETAIL
103+
104+
105+
}} // namespace boost::geometry
106+
107+
#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_COLOCATE_CLUSTERS_HPP

include/boost/geometry/algorithms/detail/overlay/get_clusters.hpp

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -39,9 +39,8 @@ struct sweep_equal_policy
3939
template <typename T>
4040
static inline T threshold()
4141
{
42-
// Tuned by cases of #1081. It should just be one epsilon to distinguish between
43-
// different turn-points and real colocated clusters.
44-
return T(1);
42+
// Points within some epsilons are considered as equal.
43+
return T(100);
4544
}
4645
public:
4746
// Returns true if point are considered equal (within an epsilon)

include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include <boost/geometry/core/point_order.hpp>
3030
#include <boost/geometry/algorithms/detail/overlay/cluster_info.hpp>
3131
#include <boost/geometry/algorithms/detail/overlay/do_reverse.hpp>
32+
#include <boost/geometry/algorithms/detail/overlay/colocate_clusters.hpp>
3233
#include <boost/geometry/algorithms/detail/overlay/get_clusters.hpp>
3334
#include <boost/geometry/algorithms/detail/overlay/get_ring.hpp>
3435
#include <boost/geometry/algorithms/detail/overlay/is_self_turn.hpp>
@@ -52,22 +53,22 @@ namespace boost { namespace geometry
5253
namespace detail { namespace overlay
5354
{
5455

56+
// Removes clusters which have only one point left, or are empty.
5557
template <typename Turns, typename Clusters>
5658
inline void remove_clusters(Turns& turns, Clusters& clusters)
5759
{
58-
typename Clusters::iterator it = clusters.begin();
60+
auto it = clusters.begin();
5961
while (it != clusters.end())
6062
{
6163
// Hold iterator and increase. We can erase cit, this keeps the
6264
// iterator valid (cf The standard associative-container erase idiom)
63-
typename Clusters::iterator current_it = it;
65+
auto current_it = it;
6466
++it;
6567

66-
std::set<signed_size_type> const& turn_indices
67-
= current_it->second.turn_indices;
68+
auto const& turn_indices = current_it->second.turn_indices;
6869
if (turn_indices.size() == 1)
6970
{
70-
signed_size_type const turn_index = *turn_indices.begin();
71+
auto const turn_index = *turn_indices.begin();
7172
turns[turn_index].cluster_id = -1;
7273
clusters.erase(current_it);
7374
}
@@ -78,18 +79,16 @@ template <typename Turns, typename Clusters>
7879
inline void cleanup_clusters(Turns& turns, Clusters& clusters)
7980
{
8081
// Removes discarded turns from clusters
81-
for (typename Clusters::iterator mit = clusters.begin();
82-
mit != clusters.end(); ++mit)
82+
for (auto& pair : clusters)
8383
{
84-
cluster_info& cinfo = mit->second;
85-
std::set<signed_size_type>& ids = cinfo.turn_indices;
86-
for (std::set<signed_size_type>::iterator sit = ids.begin();
87-
sit != ids.end(); /* no increment */)
84+
auto& cinfo = pair.second;
85+
auto& ids = cinfo.turn_indices;
86+
for (auto sit = ids.begin(); sit != ids.end(); /* no increment */)
8887
{
89-
std::set<signed_size_type>::iterator current_it = sit;
88+
auto current_it = sit;
9089
++sit;
9190

92-
signed_size_type const turn_index = *current_it;
91+
auto const turn_index = *current_it;
9392
if (turns[turn_index].discarded)
9493
{
9594
ids.erase(current_it);
@@ -98,6 +97,7 @@ inline void cleanup_clusters(Turns& turns, Clusters& clusters)
9897
}
9998

10099
remove_clusters(turns, clusters);
100+
colocate_clusters(clusters, turns);
101101
}
102102

103103
template <typename Turn, typename IdSet>

include/boost/geometry/algorithms/detail/overlay/traversal.hpp

Lines changed: 30 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -968,31 +968,43 @@ public :
968968
{
969969
turn_type const& current_turn = m_turns[turn_index];
970970

971+
bool const back_at_start_cluster
972+
= has_points
973+
&& current_turn.is_clustered()
974+
&& m_turns[start_turn_index].cluster_id == current_turn.cluster_id;
971975
if (BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection))
972976
{
973-
if (has_points)
974-
{
975-
bool const back_at_start_cluster
976-
= current_turn.is_clustered()
977-
&& m_turns[start_turn_index].cluster_id == current_turn.cluster_id;
977+
// Intersection or difference
978978

979-
if (turn_index == start_turn_index || back_at_start_cluster)
980-
{
981-
// Intersection can always be finished if returning
982-
turn_index = start_turn_index;
983-
op_index = start_op_index;
984-
return true;
985-
}
979+
if (has_points && (turn_index == start_turn_index || back_at_start_cluster))
980+
{
981+
// Intersection can always be finished if returning
982+
turn_index = start_turn_index;
983+
op_index = start_op_index;
984+
return true;
986985
}
987986

988987
if (! current_turn.is_clustered()
989-
&& current_turn.both(operation_intersection))
990-
{
991-
if (analyze_ii_intersection(turn_index, op_index,
988+
&& current_turn.both(operation_intersection)
989+
&& analyze_ii_intersection(turn_index, op_index,
992990
current_turn, previous_seg_id))
993-
{
994-
return true;
995-
}
991+
{
992+
return true;
993+
}
994+
}
995+
else if (turn_index == start_turn_index || back_at_start_cluster)
996+
{
997+
// Union or buffer: cannot return immediately to starting turn, because it then
998+
// might miss a formed multi polygon with a touching point.
999+
auto const& current_op = current_turn.operations[op_index];
1000+
signed_size_type const next_turn_index = current_op.enriched.get_next_turn_index();
1001+
bool const to_other_turn = next_turn_index >= 0 && m_turns[next_turn_index].cluster_id != current_turn.cluster_id;
1002+
if (! to_other_turn)
1003+
{
1004+
// Return to starting point
1005+
turn_index = start_turn_index;
1006+
op_index = start_op_index;
1007+
return true;
9961008
}
9971009
}
9981010

test/algorithms/buffer/buffer_multi_polygon.cpp

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -540,7 +540,10 @@ void test_all()
540540

541541
test_one<multi_polygon_type, polygon_type>("rt_p11", rt_p11, join_miter, end_flat, 28.7426, 1.0);
542542
test_one<multi_polygon_type, polygon_type>("rt_p12", rt_p12, join_miter, end_flat, 22.5711, 1.0);
543+
544+
// Needs centroid of cluster turn points
543545
test_one<multi_polygon_type, polygon_type>("rt_p13", rt_p13, join_miter, end_flat, 19.9142, 1.0);
546+
544547
test_one<multi_polygon_type, polygon_type>("rt_p14", rt_p14, join_miter, end_flat, 20.8284, 1.0);
545548
test_one<multi_polygon_type, polygon_type>("rt_p15", rt_p15, join_miter, end_flat, 23.6569, 1.0);
546549
test_one<multi_polygon_type, polygon_type>("rt_p16", rt_p16, join_miter, end_flat, 23.4853, 1.0);
@@ -631,11 +634,9 @@ void test_all()
631634
test_one<multi_polygon_type, polygon_type>("nores_6061", nores_6061, join_round32, end_flat, 39.7371, 1.0);
632635
test_one<multi_polygon_type, polygon_type>("nores_37f6", nores_37f6, join_round32, end_flat, 26.5339, 1.0);
633636

634-
#if defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES)
635-
// Fails since get_cluster works with an epsilon of 1 (instead of 1000 before).
636-
// With 3 it still succeeds but that causes regression in issue #issue_1081b
637+
// Needs an espilon in get_cluster of 3 or higher
637638
test_one<multi_polygon_type, polygon_type>("nores_1ea1", nores_1ea1, join_round32, end_flat, 28.9755, 1.0);
638-
#endif
639+
639640
test_one<multi_polygon_type, polygon_type>("nores_804e", nores_804e, join_round32, end_flat, 26.4503, 1.0);
640641
test_one<multi_polygon_type, polygon_type>("nores_51c6", nores_51c6, join_round32, end_flat, 20.2419, 1.0);
641642
test_one<multi_polygon_type, polygon_type>("nores_e5f3", nores_e5f3, join_round32, end_flat, 14.5503, 1.0);
@@ -651,7 +652,7 @@ void test_all()
651652
#if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES)
652653
// Erroneous case with rescaling
653654
test_one<multi_polygon_type, polygon_type>("res_8618", res_8618, join_round32, end_flat, 48.1085, 1.0);
654-
#endif
655+
#endif
655656
test_one<multi_polygon_type, polygon_type>("res_3b4d", res_3b4d, join_round32, end_flat, 48.4739, 1.0);
656657

657658
test_one<multi_polygon_type, polygon_type>("neighbouring_small",

test/algorithms/overlay/get_clusters.cpp

Lines changed: 5 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -101,13 +101,11 @@ int test_main(int, char* [])
101101
test_get_clusters<dp>();
102102
test_get_clusters<ep>();
103103

104-
// This constant relates to the threshold in get_clusters,
105-
// which is now return T(1) (earlier it was 1000)
106-
double const multiplier = 1.0 / 1000.0;
107-
108-
test_get_clusters_border_cases<fp>(1.0e-4 * multiplier);
109-
test_get_clusters_border_cases<dp>(1.0e-13 * multiplier);
110-
test_get_clusters_border_cases<ep>(1.0e-16 * multiplier);
104+
// These constant relate to the threshold in get_clusters.hpp,
105+
// and the used floating point type.
106+
test_get_clusters_border_cases<fp>(1.0e-5);
107+
test_get_clusters_border_cases<dp>(1.0e-14);
108+
test_get_clusters_border_cases<ep>(1.0e-17);
111109

112110
return 0;
113111
}

test/algorithms/overlay/multi_overlay_cases.hpp

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1574,6 +1574,12 @@ static std::string issue_930[2] =
15741574
"MULTIPOLYGON(((-10 2,5 3,20 2,-10 2)))"
15751575
};
15761576

1577+
static std::string issue_1109[2] =
1578+
{
1579+
"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)))",
1580+
"MULTIPOLYGON(((0 -88,0 -115.40000152587890625,-10 -88,0 -88)))"
1581+
};
1582+
15771583
static std::string bug_21155501[2] =
15781584
{
15791585
"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)))",

test/algorithms/overlay/overlay_cases.hpp

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1101,6 +1101,12 @@ static std::string issue_1081c[2] =
11011101
"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))"
11021102
};
11031103

1104+
static std::string issue_1108[2] =
1105+
{
1106+
"POLYGON((17 -2,18 -1.269679093926235014,12 0,22 0,17 -2))",
1107+
"POLYGON((22 1,22 0,14 0,18 -1.2696790939262529996,12 0,22 1))"
1108+
};
1109+
11041110
static std::string ggl_list_20120229_volker[3] =
11051111
{
11061112
"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))",

test/algorithms/set_operations/union/union.cpp

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -416,7 +416,10 @@ void test_areal()
416416
TEST_UNION(ticket_10108_a, count_set(1, 2), 0, 8, 0.0435229);
417417
TEST_UNION(ticket_10108_b, count_set(1, 2), 0, 10, 2424.3449);
418418

419+
#if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES)
420+
// With rescaling, there is a dependency on cluster tolerance, which alters the result.
419421
TEST_UNION(ticket_10866, 1, 0, 14, 332760303.5);
422+
#endif
420423

421424
TEST_UNION(ticket_11725, 1, 1, 10, 7.5);
422425

@@ -456,6 +459,9 @@ void test_areal()
456459
TEST_UNION(issue_1081c, 1, 1, -1, 2338.08);
457460
TEST_UNION_REV(issue_1081c, 1, 1, -1, 2338.08);
458461

462+
TEST_UNION(issue_1108, 1, 0, -1, 12.1742);
463+
TEST_UNION_REV(issue_1108, 1, 0, -1, 12.1742);
464+
459465
{
460466
// Rescaling produces an invalid result
461467
ut_settings settings;

test/algorithms/set_operations/union/union_multi.cpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -445,6 +445,8 @@ void test_areal()
445445
TEST_UNION(issue_888_34, 15, 0, -1, 0.3017459);
446446
TEST_UNION(issue_888_37, 52, 3, -1, 0.4033294);
447447

448+
TEST_UNION(issue_1109, 2, 0, -1, 3946.5);
449+
448450
// One or two polygons, the ideal case is 1
449451
TEST_UNION(mail_2019_01_21_johan, count_set(1, 2), 0, -1, 0.00058896);
450452

0 commit comments

Comments
 (0)