forked from boostorg/geometry
-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_clusters.hpp
220 lines (188 loc) · 6.49 KB
/
get_clusters.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2021 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_GET_CLUSTERS_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_GET_CLUSTERS_HPP
#include <algorithm>
#include <map>
#include <vector>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/algorithms/detail/overlay/approximately_equals.hpp>
#include <boost/geometry/algorithms/detail/overlay/cluster_info.hpp>
#include <boost/geometry/algorithms/detail/overlay/get_ring.hpp>
#include <boost/geometry/algorithms/detail/recalculate.hpp>
#include <boost/geometry/policies/robustness/rescale_policy_tags.hpp>
#include <boost/range/value_type.hpp>
#include <boost/geometry/util/math.hpp>
#define BOOST_GEOMETRY_USE_RESCALING_IN_GET_CLUSTERS
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DETAIL
namespace detail { namespace overlay
{
template <typename Tag = no_rescale_policy_tag, bool Integral = false>
struct sweep_equal_policy
{
public:
// Returns true if point are considered equal (within an epsilon)
template <typename P>
static inline bool equals(P const& p1, P const& p2)
{
using coor_t = typename coordinate_type<P>::type;
return approximately_equals(p1, p2, common_approximately_equals_epsilon<coor_t>::value());
}
template <typename T>
static inline bool exceeds(T value)
{
T const limit = T(1) / common_approximately_equals_epsilon<T>::value();
return value > limit;
}
};
template <>
struct sweep_equal_policy<no_rescale_policy_tag, true>
{
template <typename P>
static inline bool equals(P const& p1, P const& p2)
{
return geometry::get<0>(p1) == geometry::get<0>(p2)
&& geometry::get<1>(p1) == geometry::get<1>(p2);
}
template <typename T>
static inline bool exceeds(T value)
{
return value > 0;
}
};
#ifdef BOOST_GEOMETRY_USE_RESCALING_IN_GET_CLUSTERS
template <>
struct sweep_equal_policy<rescale_policy_tag, true>
{
template <typename P>
static inline bool equals(P const& p1, P const& p2)
{
// Neighbouring cells in the "integer grid" are considered as equal
return math::abs(geometry::get<0>(p1) - geometry::get<0>(p2)) <= 1
&& math::abs(geometry::get<1>(p1) - geometry::get<1>(p2)) <= 1;
}
template <typename T>
static inline bool exceeds(T value)
{
return value > 1;
}
};
#endif
template <typename Point>
struct turn_with_point
{
std::size_t turn_index;
Point pnt;
};
template <typename Cluster, typename Point>
struct cluster_with_point
{
Cluster cluster;
Point pnt;
};
// Use a sweep algorithm to detect clusters
template
<
typename Turns,
typename Clusters,
typename RobustPolicy
>
inline void get_clusters(Turns& turns, Clusters& clusters,
RobustPolicy const& robust_policy)
{
using turn_type = typename boost::range_value<Turns>::type;
using cluster_type = typename Clusters::mapped_type;
#ifdef BOOST_GEOMETRY_USE_RESCALING_IN_GET_CLUSTERS
// For now still use robust points for rescaled, otherwise points do not match
using point_type = typename geometry::robust_point_type
<
typename turn_type::point_type,
RobustPolicy
>::type;
#else
using point_type = typename turn_type::point_type;
#endif
sweep_equal_policy
<
typename rescale_policy_type<RobustPolicy>::type,
std::is_integral<typename coordinate_type<point_type>::type>::value
> equal_policy;
std::vector<turn_with_point<point_type>> points;
std::size_t turn_index = 0;
for (auto const& turn : turns)
{
if (! turn.discarded)
{
#ifdef BOOST_GEOMETRY_USE_RESCALING_IN_GET_CLUSTERS
point_type pnt;
geometry::recalculate(pnt, turn.point, robust_policy);
points.push_back({turn_index, pnt});
#else
points.push_back({turn_index, turn.point});
#endif
}
turn_index++;
}
// Sort the points from top to bottom
std::sort(points.begin(), points.end(), [](auto const& e1, auto const& e2)
{
return geometry::get<1>(e1.pnt) > geometry::get<1>(e2.pnt);
});
// The output vector will be sorted from bottom too
std::vector<cluster_with_point<cluster_type, point_type>> clustered_points;
// Compare points with each other. Performance is O(n log(n)) because of the sorting.
for (auto it1 = points.begin(); it1 != points.end(); ++it1)
{
// Inner loop, iterates until it exceeds coordinates in y-direction
for (auto it2 = it1 + 1; it2 != points.end(); ++it2)
{
auto const d = geometry::get<1>(it1->pnt) - geometry::get<1>(it2->pnt);
if (equal_policy.exceeds(d))
{
// Points at this y-coordinate or below cannot be equal
break;
}
if (equal_policy.equals(it1->pnt, it2->pnt))
{
std::size_t cindex = 0;
// Most recent clusters (with this y-value) are at the bottom
// therefore we can stop as soon as the y-value is out of reach (TODO)
bool found = false;
for (auto cit = clustered_points.begin();
cit != clustered_points.end(); ++cit)
{
found = equal_policy.equals(cit->pnt, it1->pnt);
if (found)
{
break;
}
cindex++;
}
// Add new cluster
if (! found)
{
cindex = clustered_points.size();
cluster_type newcluster;
clustered_points.push_back({newcluster, it1->pnt});
}
clustered_points[cindex].cluster.turn_indices.insert(it1->turn_index);
clustered_points[cindex].cluster.turn_indices.insert(it2->turn_index);
}
}
}
// Convert to map
signed_size_type cluster_id = 1;
for (auto& trace : clustered_points)
{
clusters[cluster_id++] = trace.cluster;
}
}
}} // namespace detail::overlay
#endif //DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_GET_CLUSTERS_HPP