overlay/graph/select_edge.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2025 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // Use, modification and distribution is subject to the Boost Software License,
0006 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0007 // http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SELECT_EDGE_HPP
0010 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SELECT_EDGE_HPP
0011
0012 #include <boost/core/ignore_unused.hpp>
0013 #include <boost/geometry/algorithms/detail/overlay/approximately_equals.hpp>
0014 #include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp>
0015 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
0016 #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp>
0017 #include <boost/geometry/algorithms/detail/overlay/turn_operation_id.hpp>
0018 #include <boost/geometry/algorithms/detail/overlay/graph/node_util.hpp>
0019 #include <boost/geometry/algorithms/detail/overlay/graph/select_toi_by_incoming.hpp>
0020
0021 #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSE_GRAPH)
0022 #include <boost/geometry/io/wkt/wkt.hpp>
0023 #endif
0024
0025 namespace boost { namespace geometry
0026 {
0027
0028 #ifndef DOXYGEN_NO_DETAIL
0029 namespace detail { namespace overlay
0030 {
0031
0032 template <typename Point>
0033 struct edge_and_side
0034 {
0035     turn_operation_id toi{0};
0036     Point point{};
0037     int side{0};
0038 };
0039
0040 template
0041 <
0042     bool Reverse1,
0043     bool Reverse2,
0044     overlay_type OverlayType,
0045     typename Geometry1,
0046     typename Geometry2,
0047     typename Turns,
0048     typename Clusters,
0049     typename Strategy
0050 >
0051 struct edge_selector
0052 {
0053 private:
0054     static constexpr operation_type target_operation = operation_from_overlay<OverlayType>::value;
0055     using point_type = typename Turns::value_type::point_type;
0056     using edge_type = edge_and_side<point_type>;
0057     using edges_type = std::vector<edge_type>;
0058
0059     // Use the coordinate type, but if it is too small (e.g. std::int16), use a double
0060     using coor_type = typename geometry::select_most_precise
0061         <
0062             geometry::coordinate_type_t<point_type>,
0063             double
0064         >::type;
0065
0066     // Walks over a ring to get the point after the turn.
0067     // The turn can be located at the very end of a segment.
0068     // Therefore it can be the first point on the next segment.
0069     template <typename Operation>
0070     point_type walk_to_point_after_turn(Operation const& op, point_type const& turn_point) const
0071     {
0072         static const coor_type tolerance
0073             = common_approximately_equals_epsilon_multiplier<coor_type>::value();
0074         int offset = 1;
0075         point_type point;
0076         do
0077         {
0078             geometry::copy_segment_point<Reverse1, Reverse2>(m_geometry1, m_geometry2,
0079                 op.seg_id, offset, point);
0080             ++offset;
0081         } while (approximately_equals(point, turn_point, tolerance) && offset < 10);
0082         return point;
0083     }
0084
0085     // Compares and returns true for the left most operation.
0086     // p1 is the point before the current turn.
0087     // p2 is the current turn.
0088     // So (p1, p2) together define the direction of the segment.
0089     bool select_collinear_target_edge(edge_type const& a, edge_type const& b) const
0090     {
0091         auto const& turn_a = m_turns[a.toi.turn_index];
0092         auto const& turn_b = m_turns[b.toi.turn_index];
0093         auto const& op_a = turn_a.operations[a.toi.operation_index];
0094         auto const& op_b = turn_b.operations[b.toi.operation_index];
0095
0096         auto const target_a = get_node_id(m_turns, op_a.enriched.travels_to_ip_index);
0097         auto const target_b = get_node_id(m_turns, op_b.enriched.travels_to_ip_index);
0098
0099         auto const& other_op_a = turn_a.operations[1 - a.toi.operation_index];
0100         auto const& other_op_b = turn_b.operations[1 - b.toi.operation_index];
0101
0102         if (other_op_a.enriched.travels_to_ip_index == -1)
0103         {
0104             return true;
0105         }
0106         if (other_op_b.enriched.travels_to_ip_index == -1)
0107         {
0108             return false;
0109         }
0110
0111         auto const other_target_a = get_node_id(m_turns, other_op_a.enriched.travels_to_ip_index);
0112         auto const other_target_b = get_node_id(m_turns, other_op_b.enriched.travels_to_ip_index);
0113
0114         if (target_b == other_target_a || target_b == other_target_b)
0115         {
0116             // The second edge goes via one of the targets of the first
0117             return false;
0118         }
0119         if (target_a == other_target_a || target_a == other_target_b)
0120         {
0121             // Vice versa
0122             return true;
0123         }
0124
0125         return true;
0126     }
0127
0128     void report(const char* caption, edges_type const& edges,
0129         point_type const& p1, point_type const& p2) const
0130     {
0131 #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSE_GRAPH)
0132         std::cout << " *** Sorted edges " << caption
0133         << " from " << geometry::wkt(p1) << " to " << geometry::wkt(p2)
0134         << std::endl;
0135         for (auto const& item : edges)
0136         {
0137             auto const& op = m_turns[item.toi.turn_index].operations[item.toi.operation_index];
0138             std::cout << "  -> " << item.toi
0139                 << " to " << op.enriched.travels_to_ip_index
0140                 << " side: " << item.side
0141                 << std::endl;
0142         }
0143 #endif
0144     }
0145
0146     turn_operation_id select_by_side(edges_type& edges, point_type const& p1, point_type const& p2) const
0147     {
0148         // Select point and calculate side for each edge
0149         auto const side_strategy = m_intersection_strategy.side();
0150         for (auto& edge : edges)
0151         {
0152             auto const& op = m_turns[edge.toi.turn_index].operations[edge.toi.operation_index];
0153             edge.point = walk_to_point_after_turn(op, p2);
0154             edge.side = side_strategy.apply(p1, p2, edge.point);
0155         }
0156
0157         // Sort by side (with respect to segment [p1..p2]) (TEMPORARY: and then by toi)
0158         // Right = -1 will come first. Left = 1 will come last.
0159         // This works for both union and intersection operations, because it should always
0160         // take the right turn (even in uu in buffer/union).
0161         std::sort(edges.begin(), edges.end(), [](auto const& a, auto const& b)
0162         {
0163             return std::tie(a.side, a.toi) < std::tie(b.side, b.toi);
0164         });
0165
0166         report("by side", edges, p1, p2);
0167
0168         if (edges.size() == 1 || (edges.size() > 1 && edges.front().side != edges[1].side))
0169         {
0170             return edges.front().toi;
0171         }
0172
0173         if (edges.front().side != edges.back().side)
0174         {
0175             // Remove all edges with different side than the first
0176             auto it = std::find_if(edges.begin() + 1, edges.end(), [&](auto const& item)
0177             {
0178                 return item.side != edges.front().side;
0179             });
0180             edges.erase(it, edges.end());
0181         }
0182
0183         if (edges.front().side == 0)
0184         {
0185             // Select for collinearity (it makes no sense to sort on mutual side)
0186             auto compare = [&](edge_type const& a, edge_type const& b) -> bool
0187             {
0188                 return select_collinear_target_edge(a, b);
0189             };
0190             std::sort(edges.begin(), edges.end(), compare);
0191             return edges.front().toi;
0192         }
0193
0194         // Phase 2, sort by mutual side, of the edges having the front edge's side.
0195         auto compare_one_side = [&](auto const& a, auto const& b) -> bool
0196         {
0197             // Calculating one side is enough. Either both are 0, or they are opposite.
0198             int const side = side_strategy.apply(p2, a.point, b.point);
0199             return side == 1;
0200         };
0201
0202         std::sort(edges.begin(), edges.end(), compare_one_side);
0203
0204         report("by mutual side", edges, p1, p2);
0205
0206         return edges.front().toi;
0207     }
0208
0209 public:
0210
0211     edge_selector(Geometry1 const& m_geometry1, Geometry2 const& m_geometry2,
0212         Turns const& m_turns, Clusters const& clusters,
0213         Strategy const& strategy)
0214         : m_geometry1(m_geometry1)
0215         , m_geometry2(m_geometry2)
0216         , m_turns(m_turns)
0217         , m_clusters(clusters)
0218         , m_intersection_strategy(strategy)
0219     {}
0220
0221     // Select one operation which is the leftmost or rightmost operation.
0222     // p1 is the point before the current turn.
0223     // p2 is the current turn.
0224     // So (p1, p2) together define the direction of the segment.
0225     turn_operation_id select_target_edge(set_of_tois const& turn_operation_ids,
0226             point_type const& p1, point_type const& p2) const
0227     {
0228         if (turn_operation_ids.empty())
0229         {
0230             return {};
0231         }
0232         if (turn_operation_ids.size() == 1)
0233         {
0234             return *turn_operation_ids.begin();
0235         }
0236
0237         edges_type edges;
0238         edges.reserve(turn_operation_ids.size());
0239         for (auto const& toi : turn_operation_ids)
0240         {
0241             edges.emplace_back(edge_type{toi});
0242         }
0243
0244         // Verification function for clusters: if it is clustered, all should come from one cluster.
0245         auto assert_one_cluster = [&]() -> bool
0246         {
0247             auto const& turn0 = m_turns[edges[0].toi.turn_index];
0248             auto const cluster_id = turn0.cluster_id;
0249             for (auto const& toi : turn_operation_ids)
0250             {
0251                 auto const& turn = m_turns[toi.turn_index];
0252                 if (turn.cluster_id != cluster_id)
0253                 {
0254                     return false;
0255                 }
0256             }
0257             return true;
0258         };
0259
0260         boost::ignore_unused(assert_one_cluster);
0261
0262         // It often happens there are just two collinear edges.
0263         // If they travel to the same target, take either.
0264         if (edges.size() == 2)
0265         {
0266             auto const& turn0 = m_turns[edges[0].toi.turn_index];
0267             auto const& turn1 = m_turns[edges[1].toi.turn_index];
0268             auto const& op0 = turn0.operations[edges[0].toi.operation_index];
0269             auto const& op1 = turn1.operations[edges[1].toi.operation_index];
0270             if (op0.operation == operation_continue
0271                 && op1.operation == operation_continue
0272                 && op0.enriched.travels_to_ip_index == op1.enriched.travels_to_ip_index)
0273             {
0274                 return edges.front().toi;
0275             }
0276
0277             if (target_operation == operation_union
0278                 && turn0.is_clustered()
0279                 && op0.operation == operation_union
0280                 && op1.operation == operation_union
0281                 && op0.enriched.rank == op1.enriched.rank)
0282             {
0283                 // Because it is clustered, and all operations come from the same cluster,
0284                 // the rank can be used, which is more efficient.
0285                 BOOST_GEOMETRY_ASSERT(assert_one_cluster());
0286
0287                 turn_operation_id result;
0288                 if (select_toi_for_union(result, op0, op1, edges[0].toi, edges[1].toi, m_turns))
0289                 {
0290                     return result;
0291                 }
0292
0293                 bool const better = is_better_collinear_for_union(
0294                         op0, op1, edges.front().toi, edges.back().toi);
0295                 return better ? edges.front().toi : edges.back().toi;
0296             }
0297         }
0298
0299         return select_by_side(edges, p1, p2);
0300     }
0301
0302 private:
0303     Geometry1 const& m_geometry1;
0304     Geometry2 const& m_geometry2;
0305     Turns const& m_turns;
0306     Clusters const& m_clusters;
0307     Strategy const& m_intersection_strategy;
0308 };
0309
0310 }} // namespace detail::overlay
0311 #endif // DOXYGEN_NO_DETAIL
0312
0313 }} // namespace boost::geometry
0314
0315 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SELECT_EDGE_HPP