From 43146530b0921d8276de0ae04c57cd4dc86a8bfc Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Mon, 2 Mar 2026 21:48:51 +0100 Subject: [PATCH] arrange polygons: Alternative (unused) perimiter approach; simpler topology handling; projection-based clean-up --- src/svgfill/src/arrange_polygons.cpp | 545 ++++++++++++++++++++++++++- 1 file changed, 536 insertions(+), 9 deletions(-) diff --git a/src/svgfill/src/arrange_polygons.cpp b/src/svgfill/src/arrange_polygons.cpp index 0e77a28886..4fa20c68dc 100644 --- a/src/svgfill/src/arrange_polygons.cpp +++ b/src/svgfill/src/arrange_polygons.cpp @@ -1,4 +1,4 @@ -// #define SVGFILL_DEBUG +#define SVGFILL_DEBUG // #define SVGFILL_MAIN #ifndef SVGFILL_MAIN @@ -38,6 +38,7 @@ typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Polygon_2 Polygon_2; typedef CGAL::Polygon_with_holes_2 Polygon_with_holes_2; typedef K::Point_2 Point_2; +typedef K::Vector_2 Vector_2; typedef K::Segment_2 Segment_2; typedef std::vector Polygon_list; typedef CGAL::Arr_segment_traits_2 Traits_2; @@ -694,6 +695,50 @@ class SegmentLookup { return std::make_pair(input_it, closest); }; + std::vector> n_closest_input_segments(const Segment_2& e, size_t n = 2) const { + auto mid = CGAL::ORIGIN + ((e.source() - CGAL::ORIGIN) + (e.target() - CGAL::ORIGIN)) / 2; + + std::vector>::iterator> cands; + cands.reserve(64); + + auto mid3 = CGAL::Point_3(mid.x(), mid.y(), 0); + + for (int i = -1; i <= 3; ++i) { + cands.clear(); + double r = std::pow(10.0, i); + auto midbb = mid3.bbox(); + CGAL::Bbox_3 box(midbb.xmin() - r, midbb.ymin() - r, -1.0, midbb.xmax() + r, midbb.ymax() + r, +1.0); + tree_.all_intersected_primitives(box, std::back_inserter(cands)); + if (cands.size() >= n) { + break; + } + } + + if (cands.empty()) { + return {}; + } + + std::vector>> scored; + scored.reserve(cands.size()); + for (auto it : cands) { + const auto& s3 = *it; + Segment_2 s2(Point_2(s3.source().x(), s3.source().y()), Point_2(s3.target().x(), s3.target().y())); + scored.emplace_back(CGAL::squared_distance(mid, s2), s2); + } + + std::sort(scored.begin(), scored.end(), [](auto& a, auto& b) { return a.first < b.first; }); + if (scored.size() > n) { + scored.resize(n); + } + + std::vector> out; + out.reserve(scored.size()); + for (auto& p : scored) { + out.push_back(p.second); + } + return out; + } + private: using TreeTraits = CGAL::AABB_traits>::iterator>>; using Tree = CGAL::AABB_tree; @@ -934,6 +979,69 @@ void eliminate_colinear_vertices(Graph2D& G) { } } +struct Ccw_radial_sort { + Point_2 c; + explicit Ccw_radial_sort(const Point_2& center) : c(center) {} + + bool operator()(const Point_2& a, const Point_2& b) const { + const Vector_2 va = a - c; + const Vector_2 vb = b - c; + + // Only left-turn is not sufficient because we should not wrap around, + // but rather start from e.g positive x-axis and then sort CCW. + // Therefore top-half plane always comes before bottom-half plane. + const bool ua = va.y() == 0 ? va.x() > 0 : va.y() > 0; + const bool ub = vb.y() == 0 ? vb.x() > 0 : vb.y() > 0; + if (ua != ub) { + return ua; + } + + if (CGAL::collinear(c, a, b)) { + // Nearer first so that original polygon edges are likely retained + // (not sure if it matters). + return va.squared_length() < vb.squared_length(); + } + + // This is a less functor, so we return true if c,a,b is a left turn, which means that a is CCW before b + return CGAL::left_turn(c, a, b); + } +}; + +void build_radial_neighbour_map(const std::vector& polygons, double radius, std::map>& neighbour_map) { + for (auto& poly : polygons) { + for (auto it = poly.edges_begin(); it != poly.edges_end(); ++it) { + auto source = it->source(); + auto target = it->target(); + neighbour_map[source].push_back(target); + neighbour_map[target].push_back(source); + } + } + + // Box_intersection_d package to find close vertices and connect them as well + typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + std::vector boxes; + for (auto& poly : polygons) { + for (auto it = poly.vertices_begin(); it != poly.vertices_end(); ++it) { + const auto pb = it->bbox(); + boxes.emplace_back( + CGAL::Bbox_2(pb.xmin() - radius, pb.ymin() - radius, pb.xmax() + radius, pb.ymax() + radius), + *it); + } + } + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), [&](const Box& a, const Box& b) { + if ((a.handle() - b.handle()).squared_length() <= (radius * radius)) { + neighbour_map[a.handle()].push_back(b.handle()); + neighbour_map[b.handle()].push_back(a.handle()); + } + }); + + // radial sort + for (auto& p : neighbour_map) { + auto& nb = p.second; + std::sort(nb.begin(), nb.end(), Ccw_radial_sort(p.first)); + } +} + void edge_slide(Graph2D& G) { std::list> edges_to_remove, edges_to_insert; @@ -1043,6 +1151,7 @@ std::list> extend_end_vertices_based_on_input( if (segment_to_input_facet.find(*q)->second.size() == 2) { for (auto& bnd : inner_offset) { // if point M is contained in bnd interior: + // if (!bnd.has_on_unbounded_side(M)) { if (bnd.has_on_bounded_side(M)) { auto& incoming = *it->second.begin(); // create ray incoming -> M @@ -1067,6 +1176,9 @@ std::list> extend_end_vertices_based_on_input( } if (closest_intersection_point) { + constructed_segments.push_front({M, *closest_intersection_point}); + break; +#if 0 Graph2D GGG(bnd); GGG.refine(*GGG.query(*closest_intersection_point, 0.01), *closest_intersection_point); @@ -1104,11 +1216,15 @@ std::list> extend_end_vertices_based_on_input( handled_as_graph_path = true; break; } +#endif + } else { + std::cerr << "Warning: no intersection found when extending end vertex, this will likely result in invalid topology" << std::endl; } } } } +#if 0 if (!handled_as_graph_path) { // else we choose to map point to the midpoint of the found two close points. @@ -1150,6 +1266,7 @@ std::list> extend_end_vertices_based_on_input( constructed_segments.push_front({avg, Q}); constructed_segments.push_front({avg, R}); } +#endif } } @@ -1238,6 +1355,372 @@ void fuse_corridor_halves_with_input(Arrangement_2& arr, Graph2D& G, SegmentL } } +#include + +class Segment_2_less { + public: + bool operator()(const Segment_2& a, const Segment_2& b) const { + if (a.source() != b.source()) { + return a.source() < b.source(); + } + return a.target() < b.target(); + } +}; + +void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { + using SK = CGAL::Simple_cartesian; + CGAL::Cartesian_converter C{}; + + auto other = [](const Segment_2& e, const Point_2& v) { + return (e.source() == v) ? e.target() : e.source(); + }; + + std::set edges; + for (auto he = arr.edges_begin(); he != arr.edges_end(); ++he) { + auto a = he->source()->point(); + auto b = he->target()->point(); + if (a < b) { + edges.insert({a, b}); + } else { + edges.insert({b, a}); + } + } + + auto edge_badness = [&](const Segment_2& e) -> double { + auto closest = segment_lookup.n_closest_input_segments(e, 2); + if (closest.size() != 2) { + throw std::runtime_error("Unable to locate two nearby edges"); + } + + auto get_dir = [&](const Segment_2& s) { + auto a = C(s.source()); + auto b = C(s.target()); + SK::Vector_2 v = b - a; + double l = std::sqrt(v.squared_length()); + if (l <= 1e-12) { + return std::make_pair(SK::Vector_2(0, 0), 0.); + } + return std::make_pair(v / l, l); + }; + + auto [own_dir, own_length] = get_dir(e); + + auto angle = [&](const SK::Vector_2& ov) { + double d = std::abs(own_dir * ov); + if (d > 1.0) { + d = 1.0; + } + return std::acos(d); + }; + + double best = std::numeric_limits::infinity(); + for (auto& s : closest) { + auto [dv, dl] = get_dir(s); + best = std::min(best, angle(dv)); + } + return (best + 0.1) / own_length; + }; + + std::map badnesses; + for (auto& e : edges) { + badnesses[e] = edge_badness(e); + } + + double thr; + { + std::vector tmp; + tmp.reserve(badnesses.size()); + for (auto& p : badnesses) { + tmp.push_back(p.second); + } + std::nth_element(tmp.begin(), tmp.begin() + tmp.size() / 2, tmp.end()); + double med = tmp[tmp.size() / 2]; + thr = 10.0 * med; + } + + std::set bad_edges; + for (auto& p : badnesses) { + if (p.second > thr) { + bad_edges.insert(p.first); + } + } + + std::map> topo, bad_topo; + + for (auto& e : edges) { + topo[e.source()].push_back(e); + topo[e.target()].push_back(e); + } + for (auto& e : bad_edges) { + bad_topo[e.source()].push_back(e); + bad_topo[e.target()].push_back(e); + } + + std::set break_vertices; + for (auto& p : bad_topo) { + const auto& v = p.first; + if (!(topo[v].size() == 2 && bad_topo[v].size() == 2)) { + break_vertices.insert(v); + } + } + + std::set seen; + std::vector> bad_paths; + + for (auto& s : break_vertices) { + for (auto& e0 : bad_topo[s]) { + if (seen.count(e0)) { + continue; + } + + Point_2 v = s; + std::vector path; + path.push_back(s); + + auto e = e0; + while (true) { + seen.insert(e); + v = other(e, v); + path.push_back(v); + + if (break_vertices.count(v)) { + break; + } + + auto& inc = bad_topo[v]; + std::vector nxt; + nxt.reserve(2); + for (auto& ee : inc) { + if (ee != e && !seen.count(ee)) { + nxt.push_back(ee); + } + } + if (nxt.size() != 1) { + break; + } + e = nxt[0]; + } + + if (path.size() > 1) { + bad_paths.push_back(std::move(path)); + } + } + } + + std::set> to_remove; + std::vector> to_insert; + + auto dirs_from = [&](const Point_2& v, const std::set& path_edges) { + std::vector ds; + for (auto& ee : topo[v]) { + if (path_edges.count(ee) || bad_edges.count(ee)) { + continue; + } + Point_2 u = other(ee, v); + ds.push_back(v - u); + } + return ds; + }; + + auto collapse_path = [&](const std::vector& path) -> std::optional { + std::set path_edges; + for (size_t i = 0; i + 1 < path.size(); ++i) { + auto* a = &path[i]; + auto* b = &path[i + 1]; + if (*a < *b) { + path_edges.insert({*a, *b}); + } else { + path_edges.insert({*b, *a}); + } + } + + const auto& a0 = path.front(); + const auto& b0 = path.back(); + + auto das = dirs_from(a0, path_edges); + auto dbs = dirs_from(b0, path_edges); + if (das.empty() || dbs.empty()) { + return std::nullopt; + } + + K::FT best_ke = std::numeric_limits::infinity(); + std::optional best_x; + + for (auto& da : das) { + for (auto& db : dbs) { + CGAL::Ray_2 r1(a0, da); + CGAL::Ray_2 r2(b0, db); + + auto x = CGAL::intersection(r1, r2); + if (x) { + if (auto* xp = variant_get>(&*x)) { + // @todo does it matter that this is squared? + auto ke = CGAL::squared_distance(*xp, a0) + CGAL::squared_distance(*xp, b0); + if (ke < best_ke) { + best_ke = ke; + best_x = *xp; + } + } + } + } + } + + return best_x; + }; + + for (auto& path : bad_paths) { + auto x = collapse_path(path); + if (!x) { + // std::cerr << "Unable to collapse path, skipping" << std::endl; + continue; + } + + double orig_length = 0.; + for (size_t i = 0; i < path.size() - 1; ++i) { + auto& a = path[i]; + auto& b = path[i + 1]; + orig_length += std::sqrt(CGAL::to_double(CGAL::squared_distance(a, b))); + } + + double new_length = std::sqrt(CGAL::to_double((path.front() - *x).squared_length())) + std::sqrt(CGAL::to_double((path.back() - *x).squared_length())); + + if (new_length > orig_length * 2 || orig_length > new_length * 2) { + // std::cerr << "Collapsing path would increase length too much, skipping" << std::endl; + continue; + } + + std::cerr << "new_length: " << new_length << " orig_length: " << orig_length << std::endl; + + for (size_t i = 0; i < path.size(); ++i) { + auto& v = path[i]; + if (CGAL::squared_distance(v, *x) < 1.e-5) { + std::cerr << "Collapsing path would create near-duplicate vert to previous path, skipping" << std::endl; + continue; + } + } + + for (size_t i = 0; i < path.size() - 1; ++i) { + auto& a = path[i]; + auto& b = path[i + 1]; + if (a < b) { + to_remove.insert({a, b}); + } else { + to_remove.insert({b, a}); + } + } + auto s = path.front(); + auto t = path.back(); + if (s != *x) { + to_insert.push_back({s, *x}); + } + if (t != *x) { + to_insert.push_back({t, *x}); + } + } + + /* + using Walk_pl = CGAL::Arr_walk_along_line_point_location; + Walk_pl walk_pl(arr); + + for (auto& e : to_remove) { + // debug_output.write_segment(e->source()->point(), e->target()->point(), "arr_bad_remove"); + auto res = walk_pl.locate(e.first); + if (auto* v = boost::get(&res)) { + Arrangement_2::Halfedge_around_vertex_circulator first, curr; + first = curr = (*v)->incident_halfedges(); + size_t i = 0; + std::array pts; + std::array hes; + do { + Arrangement_2::Vertex_const_handle u = curr->source(); + hes[i] = curr; + pts[i++] = u->point(); + } while (++curr != first); + + + if ((*v)->point() != e.first) { + std::cerr << "Warning: unable to locate vertex for edge removal, skipping" << std::endl; + continue; + } + } else { + std::cerr << "Warning: unable to locate vertex for edge removal, skipping" << std::endl; + continue; + } + } + */ + + for (auto& e : to_remove) { + bool removed = false; + for (auto he = arr.edges_begin(); he != arr.edges_end(); ++he) { + auto a = he->source()->point(); + auto b = he->target()->point(); + if ((a == e.first && b == e.second) || (a == e.second && b == e.first)) { + CGAL::remove_edge(arr, he); + removed = true; + break; + } + } + if (!removed) { + std::cerr << "Warning: unable to locate edge for removal, skipping" << std::endl; + } + } + + for (auto& pq : to_insert) { + if (pq.first == pq.second) { + continue; + } + CGAL::insert(arr, Segment_2(pq.first, pq.second)); + // debug_output.write_segment(pq.first, pq.second, "arr_bad_insert"); + } + +} + +void remove_colinear_vertices(Arrangement_2& arr) { + std::set to_remove; + std::set> to_add; + while (true) { + bool removed_this_round = false; + for (auto it = arr.vertices_begin(); it != arr.vertices_end(); ++it) { + if (it->degree() == 2) { + Arrangement_2::Halfedge_around_vertex_circulator first, curr; + first = curr = it->incident_halfedges(); + size_t i = 0; + std::array pts; + std::array hes; + do { + Arrangement_2::Vertex_const_handle u = curr->source(); + hes[i] = curr; + pts[i++] = u->point(); + } while (++curr != first); + + using SK = CGAL::Simple_cartesian; + CGAL::Cartesian_converter C{}; + + auto a = C(pts[0]); + auto b = C(it->point()); + auto c = C(pts[1]); + + SK::Vector_2 ab = b - a; + SK::Vector_2 ac = c - a; + ab /= std::sqrt(ab.squared_length()); + ac /= std::sqrt(ac.squared_length()); + + double d = ab * ac; + if (std::acos(d) < 1e-12) { + CGAL::remove_edge(arr, hes[0]); + CGAL::remove_edge(arr, hes[1]); + CGAL::insert(arr, Segment_2(pts[0], pts[1])); + removed_this_round = true; + break; + }; + } + } + if (!removed_this_round) { + break; + } + } +} + class timer { public: class entry { @@ -1266,7 +1749,7 @@ class timer { }; void arrange_cgal_polygons(const std::vector& input_polygons_, std::vector& output_polygons, double polygon_offset_distance = -1.) { - static const double OVERLAP_RESOLUTION_DISTANCE = 1.e-2; + static const double OVERLAP_RESOLUTION_DISTANCE = 1.e-1; // even larger amount of inset so that outer perimeter is safely within all input polygons even when overlap resolution is applied // no, `1.e-2 + 1.e-5` creates issues with the outer perimeter, are there other tolerances in play? static const double OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT = 1.e-5; @@ -1325,13 +1808,12 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // // Inset-offset to remove tiny details that may cause enourmous spikes in offsets for (auto& r : input_polygons) { - smooth_polygon(-polygon_offset_distance / 10000., r); + smooth_polygon(polygon_offset_distance / 1000., r); } debug_output.write_polygons(input_polygons, "processed_input"); - SegmentLookup segment_lookup(input_polygons); - +#if 1 t0 = timer.start("outer perimeter"); // Find the outer perimeter using offset - union - negative offset @@ -1372,7 +1854,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v remove_close_points(fused_removed_close_points, 1.e-4); // Apply negative offset to get the outer perimeter polygon - auto inner_offset = create_and_convert_offset_polygon( + auto outer_perimiter = create_and_convert_offset_polygon( // Because polygon_offset is inexact, make sure our inset distance is slightly larger // std::nexttoward(-polygon_offset_distance, -std::numeric_limits::infinity()), @@ -1380,14 +1862,37 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v -polygon_offset_distance - OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT, fused_removed_close_points); - debug_output.write_polygons(inner_offset, "outer_perimiter"); + debug_output.write_polygons(outer_perimiter, "outer_perimiter"); +#else + std::map> neighbour_map; + build_radial_neighbour_map(input_polygons, polygon_offset_distance, neighbour_map); + + auto start_vertex = neighbour_map.rbegin()->first; + auto next_vertex = neighbour_map.rbegin()->second.front(); + + std::vector cycle = {start_vertex, next_vertex}; + while (cycle.back() != cycle.front()) { + const auto& incoming_from = *(cycle.rbegin() + 1); + const auto& nb = neighbour_map[cycle.back()]; + auto it = std::find(nb.begin(), nb.end(), incoming_from); + // cycle it -1 around nb + if (it == nb.begin()) { + it == nb.end() - 1; + } else { + --it; + } + cycle.push_back(*it); + } + std::vector outer_perimiter; + outer_perimiter.emplace_back(cycle.begin(), cycle.end()); +#endif t0.stop(); t0 = timer.start("corridor creation"); // Subtract original polygons from outer perimeter std::vector difference_result, difference_result_subdivided; - for (auto& i : inner_offset) { + for (auto& i : outer_perimiter) { std::vector working_copy; working_copy.emplace_back(i); @@ -1433,6 +1938,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v debug_output.write_polygons(triangular_polygons, "triangulated_corridor"); + SegmentLookup segment_lookup(input_polygons); + auto [line_graph, midpoint_to_segment, segment_to_input_facet] = build_line_graph(input_polygons, segment_lookup, triangular_polygons); for (auto& p : line_graph) { for (auto& q : p.second) { @@ -1474,7 +1981,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v t0 = timer.start("topology"); - auto segments = extend_end_vertices_based_on_input(G, midpoint_to_segment, segment_to_input_facet, inner_offset, segment_lookup); + auto segments = extend_end_vertices_based_on_input(G, midpoint_to_segment, segment_to_input_facet, outer_perimiter, segment_lookup); // Now plot the edges on an arrangement in order to find planar cycles // and merge the corridor-halves with their neighbouring input polygon @@ -1490,7 +1997,9 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v debug_output.write_segment(pq.first, pq.second, "extended_segments"); } +#if 0 // Write input polygons to arrangement_2 + // We no longer do this because we add the outer perimiter now, subdivided by the corridor network which is extended and intersected with the outer perimiter for (auto& poly : input_polygons) { for (size_t i = 0; i != poly.size(); ++i) { auto j = (i + 1) % poly.size(); @@ -1500,6 +2009,19 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); } } +#else + // Write outer perimeter to arrangement_2 + for (auto& p : outer_perimiter) { + for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { + auto source = it->source(); + auto target = it->target(); + if (source == target) { + continue; + } + CGAL::insert(arr, Segment_2(source, target)); + } + } +#endif // Just for the automatic numbering, create a full vector std::vector temp; @@ -1525,7 +2047,12 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // corridor network we know it needs to be joined with an input polygon. In that // case the edges need to be eliminated that correspond to original geometry. +#if 0 fuse_corridor_halves_with_input(arr, G, segment_lookup, input_polygons, debug_output); +#else + remove_colinear_vertices(arr); + clean_noisy_paths(arr, segment_lookup); +#endif t0.stop();