From cfea9728ac0317ddead7cfadb786fcbab8bb84cc Mon Sep 17 00:00:00 2001 From: Richard Fairhurst Date: Sun, 19 Nov 2023 18:05:02 +0000 Subject: [PATCH] Update to latest version of kleunen/boost_geometry_correct (#581) --- include/geometry/correct.hpp | 350 +++++++++++++++++++++++++---------- 1 file changed, 251 insertions(+), 99 deletions(-) diff --git a/include/geometry/correct.hpp b/include/geometry/correct.hpp index 7adb5d34..bbf1df25 100644 --- a/include/geometry/correct.hpp +++ b/include/geometry/correct.hpp @@ -18,6 +18,10 @@ #include #include +#include +#include +#include + namespace geometry { namespace impl { @@ -46,6 +50,13 @@ static inline void result_combine(C &result, T &&new_element) } } +template +static inline void result_combine_multiple(C &result, T &new_elements) +{ + for(auto &element: new_elements) + result_combine(result, std::move(element)); +} + struct pseudo_vertice_key { std::size_t index_1; @@ -84,6 +95,13 @@ struct pseudo_vertice { } }; +struct assign_policy { + static bool const include_no_turn = true; + static bool const include_degenerate = true; + static bool const include_opposite = true; + static bool const include_start_turn = true; +}; + template< typename point_t = boost::geometry::model::d2::point_xy, typename ring_t = boost::geometry::model::ring @@ -94,46 +112,59 @@ static inline void dissolve_find_intersections( std::set &start_keys) { if(ring.empty()) return; + + for(std::size_t i = 0; i < ring.size(); ++i) { + pseudo_vertices.emplace(pseudo_vertice_key(i, i, 0.0), ring[i]); + } - boost::geometry::index::rtree, std::size_t >, boost::geometry::index::quadratic<16>> index; + // Detect intersections and generate pseudo-vertices + // boost::geometry::strategies::relate::cartesian<> strategy; + typedef boost::geometry::detail::no_rescale_policy rescale_policy_type; +#if BOOST_VERSION < 107600 + typename boost::geometry::strategy::relate::services::default_strategy::type strategy; + typedef boost::geometry::detail::overlay::turn_info + < + point_t, boost::geometry::segment_ratio + > turn_info; +#else + boost::geometry::strategies::cartesian<> strategy; + typedef boost::geometry::detail::overlay::turn_info + < + point_t + > turn_info; +#endif - // Generate all by-pass intersections in the graph - // Generate a list of all by-pass intersections - pseudo_vertices.emplace(pseudo_vertice_key(ring.size() - 1, ring.size() - 1, 0.0), ring.back()); - for(std::size_t i = ring.size() - 1; i--; ) - { - pseudo_vertices.emplace(pseudo_vertice_key(i, i, 0.0), ring[i]); - boost::geometry::model::segment line_1(ring[i], ring[i + 1]); - - boost::geometry::index::query( - index, boost::geometry::index::intersects(line_1), - boost::make_function_output_iterator([&](std::pair< boost::geometry::model::segment, std::size_t > const &iter) { - - auto const &line_2 = iter.first; - auto j = iter.second; - - std::vector output; - boost::geometry::intersection(line_1, line_2, output); - - for(auto const &p: output) { - double scale_1 = boost::geometry::comparable_distance(p, ring[i]) / boost::geometry::comparable_distance(ring[i + 1], ring[i]); - double scale_2 = boost::geometry::comparable_distance(p, ring[j]) / boost::geometry::comparable_distance(ring[j + 1], ring[j]); - if(scale_1 < 1.0 && scale_2 < 1.0) { - pseudo_vertice_key key_j(j, i, scale_2); - pseudo_vertices.emplace(pseudo_vertice_key(i, j, scale_1, true), pseudo_vertice(p, key_j)); - pseudo_vertices.emplace(key_j, p); - start_keys.insert(key_j); - - pseudo_vertice_key key_i(i, j, scale_1); - pseudo_vertices.emplace(pseudo_vertice_key(j, i, scale_2, true), pseudo_vertice(p, key_i)); - pseudo_vertices.emplace(key_i, p); - start_keys.insert(key_i); - } - } - })); + std::vector turns; - index.insert(std::make_pair(boost::geometry::model::segment(ring[i], ring[i+1]), i)); - } + rescale_policy_type rescale_policy; + + boost::geometry::detail::self_get_turn_points::no_interrupt_policy policy; + boost::geometry::self_turns + < + assign_policy + >(ring, strategy, rescale_policy, turns, policy); + + for(auto const &turn: turns) { + auto p = turn.point; + auto i = std::min(turn.operations[0].seg_id.segment_index, turn.operations[1].seg_id.segment_index); + auto j = std::max(turn.operations[0].seg_id.segment_index, turn.operations[1].seg_id.segment_index); + + double offset_1 = boost::geometry::comparable_distance(p, ring[i]); + double offset_2 = boost::geometry::comparable_distance(p, ring[j]); + + double length = boost::geometry::comparable_distance(ring[i], ring[j]); + if ((offset_1 > 0 && offset_1 < length) || (offset_2 > 0 && offset_2 < length)) { + pseudo_vertice_key key_j(j, i, offset_2); + pseudo_vertices.emplace(pseudo_vertice_key(i, j, offset_1, true), pseudo_vertice(p, key_j)); + pseudo_vertices.emplace(key_j, p); + start_keys.insert(key_j); + + pseudo_vertice_key key_i(i, j, offset_1); + pseudo_vertices.emplace(pseudo_vertice_key(j, i, offset_2, true), pseudo_vertice(p, key_i)); + pseudo_vertices.emplace(key_i, p); + start_keys.insert(key_i); + } + } } // Remove invalid points (NaN) from ring @@ -183,43 +214,70 @@ static inline void correct_close(ring_t &ring) } +template< typename point_t = boost::geometry::model::d2::point_xy > +struct compare_point_less +{ + bool operator()(point_t const &a, point_t const &b) const { + if(a.x() < b.x()) return true; + if(a.x() > b.x()) return false; + return (a.y() < b.y()); + }; +}; + template< typename point_t = boost::geometry::model::d2::point_xy, typename ring_t = boost::geometry::model::ring > -static inline std::vector dissolve_generate_rings( +static inline std::vector> dissolve_generate_rings( std::map, compare_pseudo_vertice_key> &pseudo_vertices, - std::set &start_keys, + std::set const &all_start_keys, boost::geometry::order_selector order, double remove_spike_min_area = 0.0) { - std::vector result; + std::vector> result; // Generate all polygons by tracing all the intersections // Perform union to combine all polygons into single polygon again + auto start_keys = all_start_keys; while(!start_keys.empty()) { ring_t new_ring; - + // Store point in generated polygon auto push_point = [&new_ring](auto const &p) { - if(new_ring.empty() || boost::geometry::comparable_distance(new_ring.back(), p) > 0) + if(new_ring.empty() || boost::geometry::comparable_distance(new_ring.back(), p) > 0) { new_ring.push_back(p); + } }; - auto is_closed = [](ring_t &ring) { - if(ring.size() < 2) return false; - for(std::size_t i = 0; i < ring.size() - 1; ++i) { - if(boost::geometry::comparable_distance(ring[i], ring.back()) == 0) { - ring.erase(ring.begin(), ring.begin() + i); - return true; - } + // Store newly generated ring + auto push_ring = [&result, remove_spike_min_area](ring_t &new_ring) { + auto area = boost::geometry::area(new_ring); + if(std::abs(area) > remove_spike_min_area) { + result.push_back(std::make_pair(std::move(new_ring), area)); } + }; + auto i = pseudo_vertices.find(*start_keys.begin()); + + std::vector< std::pair > start_points; + start_points.push_back(std::make_pair(i->second.p, 0)); + + // Check if the outer or inner ring is closed + auto is_closed = [&new_ring, &start_points, &push_ring](point_t const &p) { + for(auto const &i: start_points) { + if(new_ring.size() > i.second+1 && boost::geometry::comparable_distance(i.first, p) == 0) { + if(i.second == 0) return true; + + // Copy the new inner ring + ring_t inner_ring(new_ring.begin() + i.second, new_ring.end()); + push_ring(inner_ring); + + // Remove the inner ring + new_ring.erase(new_ring.begin() + i.second, new_ring.end()); + } + } return false; }; - auto start_iter = pseudo_vertices.find(*start_keys.begin()); - auto i = start_iter; - do { auto const &key = i->first; auto const &value = i->second; @@ -227,7 +285,17 @@ static inline std::vector dissolve_generate_rings( // Store the point in output polygon push_point(value.p); - start_keys.erase(key); + // Remove the key from the starting keys list + auto compare_key = [&key](pseudo_vertice_key const &i) { + return (key.index_1 == i.index_1 && key.index_2 == i.index_2 && key.scale == i.scale && key.reroute == i.reroute); + }; + + start_keys.erase(key); + + // Store possible new inner ring starting point + if(all_start_keys.find(key) != all_start_keys.end()) + start_points.push_back(std::make_pair(value.p, new_ring.size() - 1)); + if(key.reroute) { // Follow by-pass i = pseudo_vertices.find(value.link); @@ -239,13 +307,10 @@ static inline std::vector dissolve_generate_rings( } // Repeat until back at starting point - } while(!is_closed(new_ring)); + } while(!is_closed(new_ring.back())); // Combine with already generated polygons - auto area = boost::geometry::area(new_ring); - if(std::abs(area) > remove_spike_min_area) { - result.push_back(std::move(new_ring)); - } + push_ring(new_ring); } return result; @@ -257,11 +322,11 @@ template< typename ring_t = boost::geometry::model::ring, typename multi_polygon_t = boost::geometry::model::multi_polygon > -static inline std::vector correct(ring_t const &ring, boost::geometry::order_selector order, double remove_spike_min_area = 0.0) +static inline std::vector> correct(ring_t const &ring, boost::geometry::order_selector order, double remove_spike_min_area = 0.0) { constexpr std::size_t min_nodes = 3; if(ring.size() < min_nodes) - return std::vector(); + return { }; std::map, compare_pseudo_vertice_key> pseudo_vertices; std::set start_keys; @@ -281,8 +346,9 @@ static inline std::vector correct(ring_t const &ring, boost::geometry::o dissolve_find_intersections(new_ring, pseudo_vertices, start_keys); if(start_keys.empty()) { - if(std::abs(boost::geometry::area(new_ring)) > remove_spike_min_area) - return { new_ring }; + double area = boost::geometry::area(new_ring); + if(std::abs(area) > remove_spike_min_area) + return { std::make_pair(new_ring, area) }; else return { }; } @@ -295,16 +361,59 @@ template< typename polygon_t = boost::geometry::model::polygon, typename multi_polygon_t = boost::geometry::model::multi_polygon > -struct combine_non_zero_winding +static inline void fill_normalize_polygons(std::vector> &input) { - inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) + for(auto &i: input) { + for(auto &poly: i.first) { + if(i.second < 0) { + std::reverse(poly.outer().begin(), poly.outer().end()); + } + } + } +} + +template< + typename point_t = boost::geometry::model::d2::point_xy, + typename polygon_t = boost::geometry::model::polygon, + typename multi_polygon_t = boost::geometry::model::multi_polygon + > +struct fill_non_zero_winding +{ + inline void operator()(std::vector> &input) const { - if(boost::geometry::area(poly) > 0) - result_combine(combined_outers, std::move(poly)); - else { - std::reverse(poly.outer().begin(), poly.outer().end()); - result_combine(combined_inners, std::move(poly)); + auto compare = [](std::pair const &a, std::pair const &b) { return std::abs(a.second) > std::abs(b.second); }; + std::sort(input.begin(), input.end(), compare); + + std::vector scores; + for(auto &mp: input) { + scores.push_back(mp.second > 0 ? 1 : -1); } + + fill_normalize_polygons(input); + + for(std::size_t i = 0; i < input.size(); ++i) { + for(std::size_t j = i + 1; j < input.size(); ++j) { + if(boost::geometry::covered_by(input[j].first, input[i].first)) { + scores[j] += scores[i]; + } + } + } + + multi_polygon_t combined_outers; + multi_polygon_t combined_inners; + + for(std::size_t i = 0; i < input.size(); ++i) { + if(scores[i] != 0) + result_combine_multiple(combined_outers, input[i].first); + else + result_combine_multiple(combined_inners, input[i].first); + } + + multi_polygon_t output; + boost::geometry::difference(combined_outers, combined_inners, output); + + input.resize(1); + input.front().first = std::move(output); } }; @@ -313,74 +422,92 @@ template< typename polygon_t = boost::geometry::model::polygon, typename multi_polygon_t = boost::geometry::model::multi_polygon > -struct combine_odd_even +struct fill_odd_even { - inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) + inline void operator()(std::vector> &input) const { - if(boost::geometry::area(poly) < 0) - std::reverse(poly.outer().begin(), poly.outer().end()); + auto compare = [](std::pair const &a, std::pair const &b) { return std::abs(a.second) < std::abs(b.second); }; + std::sort(input.begin(), input.end(), compare); + + fill_normalize_polygons(input); + + while(input.size() > 1) { + std::size_t divide_i = input.size() / 2 + input.size() % 2; + for(std::size_t i = 0; i < input.size() / 2; ++i) { + std::size_t index = i + divide_i; + if(index < input.size()) { + multi_polygon_t result; + boost::geometry::sym_difference(input[index].first, input[i].first, result); + input[i].first = std::move(result); + } + } - multi_polygon_t result; - boost::geometry::sym_difference(combined_outers, poly, result); - combined_outers = std::move(result); + input.resize(divide_i); + } } }; template< + typename fill_function_t, typename combine_function_t, + typename difference_function_t, typename point_t = boost::geometry::model::d2::point_xy, typename polygon_t = boost::geometry::model::polygon, typename ring_t = boost::geometry::model::ring, typename multi_polygon_t = boost::geometry::model::multi_polygon > -static inline void correct(polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, combine_function_t combine) +static inline void correct(polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, fill_function_t const &fill, combine_function_t const &combine, difference_function_t const &difference) { auto order = boost::geometry::point_order::value; auto outer_rings = correct(input.outer(), order, remove_spike_min_area); - // Calculate all outers and combine them if possible - multi_polygon_t combined_outers; - multi_polygon_t combined_inners; + // Calculate all outers + std::vector> combined_outers; - for(auto &ring: outer_rings) { + for(auto &i: outer_rings) { polygon_t poly; - poly.outer() = std::move(ring); - combine(combined_outers, combined_inners, poly); + poly.outer() = std::move(i.first); + + combined_outers.push_back(std::make_pair(multi_polygon_t(), i.second)); + combined_outers.back().first.push_back(std::move(poly)); } + // fill the collected outers and combine into single multi_polygon + fill(combined_outers); + // Calculate all inners and combine them if possible + multi_polygon_t combined_inners; for(auto const &ring: input.inners()) { polygon_t poly; poly.outer() = std::move(ring); multi_polygon_t new_inners; - correct(poly, new_inners, remove_spike_min_area, combine); - - for(auto &poly: new_inners) { - result_combine(combined_inners, std::move(poly)); - } + correct(poly, new_inners, remove_spike_min_area, fill, combine, difference); + combine(combined_inners, new_inners); } // Cut out all inners from all the outers - boost::geometry::difference(combined_outers, combined_inners, output); + if(!combined_outers.empty()) { + difference(combined_outers.front().first, combined_inners, output); + } } template< + typename fill_function_t, typename combine_function_t, + typename difference_function_t, typename point_t = boost::geometry::model::d2::point_xy, typename polygon_t = boost::geometry::model::polygon, typename ring_t = boost::geometry::model::ring, typename multi_polygon_t = boost::geometry::model::multi_polygon > -static inline void correct(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, combine_function_t combine) +static inline void correct(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, fill_function_t const &fill, combine_function_t const &combine, difference_function_t const &difference) { for(auto const &polygon: input) { multi_polygon_t new_polygons; - correct(polygon, new_polygons, remove_spike_min_area, combine); - - for(auto &new_polygon: new_polygons) - result_combine(output, std::move(new_polygon)); + correct(polygon, new_polygons, remove_spike_min_area, fill, combine, difference); + combine(output, new_polygons); } } @@ -393,7 +520,11 @@ template< > static inline void correct(polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area = 0.0) { - impl::correct(input, output, remove_spike_min_area, impl::combine_non_zero_winding()); + impl::correct(input, output, remove_spike_min_area, + impl::fill_non_zero_winding(), + impl::result_combine_multiple, + boost::geometry::difference + ); } template< @@ -403,7 +534,16 @@ template< > static inline void correct_odd_even(polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area = 0.0) { - impl::correct(input, output, remove_spike_min_area, impl::combine_odd_even()); + impl::correct(input, output, remove_spike_min_area, + impl::fill_odd_even(), + [](multi_polygon_t &a, multi_polygon_t const &b) { + multi_polygon_t result; + boost::geometry::sym_difference(a, b, result); + a = std::move(result); + }, + boost::geometry::sym_difference + ); + } @@ -415,7 +555,11 @@ template< > static inline void correct(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area = 0.0) { - impl::correct(input, output, remove_spike_min_area, impl::combine_non_zero_winding()); + impl::correct(input, output, remove_spike_min_area, + impl::fill_non_zero_winding(), + impl::result_combine_multiple, + boost::geometry::difference + ); } template< @@ -426,9 +570,17 @@ template< > static inline void correct_odd_even(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area = 0.0) { - impl::correct(input, output, remove_spike_min_area, impl::combine_odd_even()); + impl::correct(input, output, remove_spike_min_area, + impl::fill_odd_even(), + [](multi_polygon_t &a, multi_polygon_t const &b) { + multi_polygon_t result; + boost::geometry::sym_difference(a, b, result); + a = std::move(result); + }, + boost::geometry::sym_difference + ); } } -#endif +#endif \ No newline at end of file