ceph/src/boost/boost/geometry/algorithms/detail/overlay/get_clusters.hpp

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