diff --git a/src/ifcopenshell-python/ifcopenshell/draw.py b/src/ifcopenshell-python/ifcopenshell/draw.py index 5f6d761ceb..962dbbb34f 100644 --- a/src/ifcopenshell-python/ifcopenshell/draw.py +++ b/src/ifcopenshell-python/ifcopenshell/draw.py @@ -42,6 +42,7 @@ WHITE = numpy.array((1.0, 1.0, 1.0)) DO_NOTHING = lambda *args: None +ARRANGE_POLYGON_SETTINGS = W.arrange_polygon_settings() if hasattr(W, 'arrange_polygon_settings') else None @dataclass class draw_settings: @@ -536,7 +537,7 @@ def main( *(tup for i, tup in enumerate(zip(path_objects, section_polies, polies)) if has_relevant_zone(i)) ) - arranged = W.arrange_polygons(polies) + arranged = W.arrange_polygons(*filter(None, (ARRANGE_POLYGON_SETTINGS,)), polies) svg_data_3 = W.polygons_to_svg(arranged, False) dom3 = parseString(svg_data_3) svg3 = dom3.childNodes[0] diff --git a/src/ifcwrap/IfcGeomWrapper.i b/src/ifcwrap/IfcGeomWrapper.i index e155c2837e..e992e4beab 100644 --- a/src/ifcwrap/IfcGeomWrapper.i +++ b/src/ifcwrap/IfcGeomWrapper.i @@ -1166,6 +1166,7 @@ ifcopenshell::geometry::taxonomy::item::ptr try_upcast(PyObject* obj0, swig_type %ignore svgfill::line_segments_to_polygons; %ignore svgfill::svg_to_polygons; %ignore svgfill::arrange_polygons; +%ignore svgfill::abstract_arrangement; %template(svg_line_segments) std::vector>; %template(svg_groups_of_line_segments) std::vector>>; @@ -1287,9 +1288,9 @@ ifcopenshell::geometry::taxonomy::item::ptr try_upcast(PyObject* obj0, swig_type } } - std::vector arrange_polygons(const std::vector& polygons) { + std::vector arrange_polygons(svgfill::arrange_polygon_settings settings, const std::vector& polygons) { std::vector r; - if (svgfill::arrange_polygons(polygons, r)) { + if (svgfill::arrange_polygons(settings, polygons, r)) { return r; } else { throw std::runtime_error("Failed to arrange polygons"); diff --git a/src/svgfill/src/arrange_polygons.cpp b/src/svgfill/src/arrange_polygons.cpp index 4fa20c68dc..ba8785a267 100644 --- a/src/svgfill/src/arrange_polygons.cpp +++ b/src/svgfill/src/arrange_polygons.cpp @@ -350,6 +350,8 @@ find_overlaps(const std::vector& polygons) { class DebugWriter { public: + DebugWriter() : enabled_(false) {} + DebugWriter(bool enabled, const std::string& filename_prefix) : enabled_(enabled) { if (enabled_) { @@ -368,6 +370,43 @@ class DebugWriter { } } + DebugWriter(const DebugWriter&) = delete; + + DebugWriter(DebugWriter&& other) noexcept + : obj(std::move(other.obj)), vi(other.vi), svg(std::move(other.svg)), enabled_(other.enabled_), last_segment_name_(std::move(other.last_segment_name_)) + { + other.enabled_ = false; + other.vi = 1; + other.last_segment_name_.clear(); + } + + DebugWriter& operator=(const DebugWriter&) = delete; + + DebugWriter& operator=(DebugWriter&& other) noexcept { + if (this == &other) { + return *this; + } + + if (enabled_) { + svg << "\n"; + obj << std::flush; + obj.close(); + svg.close(); + } + + obj = std::move(other.obj); + svg = std::move(other.svg); + vi = other.vi; + enabled_ = other.enabled_; + last_segment_name_ = std::move(other.last_segment_name_); + + other.enabled_ = false; + other.vi = 1; + other.last_segment_name_.clear(); + + return *this; + } + void write_polygon(const Polygon_2& polygon, const std::string& name) { if (enabled_) { write_polygon_to_obj_(obj, vi, true, polygon, name); @@ -387,7 +426,7 @@ class DebugWriter { obj << "l " << vi++; obj << " " << vi++ << "\n"; - svg << ""; + svg << "\n"; obj << std::flush; } @@ -468,7 +507,7 @@ class DebugWriter { } }; -void eliminate_overlaps(double OVERLAP_RESOLUTION_DISTANCE, std::vector& polygons) { +void eliminate_overlaps(DebugWriter& debug_writer, double OVERLAP_RESOLUTION_DISTANCE, std::vector& polygons) { // solve overlaps by means of subtraction // loop over overlaps and subtract the smaller polygon from the larger one @@ -576,11 +615,37 @@ void eliminate_overlaps(double OVERLAP_RESOLUTION_DISTANCE, std::vector(25, 27); + bool success = false; if ((mp1 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE, *poly2)))) { + if (is_) { + debug_writer.write_polygon(*mp1, "mp1"); + } + smooth_polygon(OVERLAP_RESOLUTION_DISTANCE / 100., *mp1); + if (is_) { + debug_writer.write_polygon(*mp1, "mp1b"); + } if ((mp2 = subtract_retain_largest(*poly1, *mp1))) { + if (is_) { + debug_writer.write_polygon(*mp2, "mp2"); + } + smooth_polygon(OVERLAP_RESOLUTION_DISTANCE / 100., *mp2); + if (is_) { + debug_writer.write_polygon(*mp2, "mp2b"); + } if ((mp3 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE * 2, *mp2)))) { + if (is_) { + debug_writer.write_polygon(*mp3, "mp3"); + } + smooth_polygon(OVERLAP_RESOLUTION_DISTANCE / 100., *mp3); + if (is_) { + debug_writer.write_polygon(*mp3, "mp3b"); + } if ((mp4 = subtract_retain_largest(*poly2, *mp3))) { + if (is_) { + debug_writer.write_polygon(*mp4, "mp4"); + } *poly1 = *mp2; *poly2 = *mp4; success = true; @@ -777,14 +842,19 @@ Polygon_with_holes_2 subdivide_polygon(double max_distance, const Polygon_with_h std::tuple< std::map>, std::map>, - std::map, std::vector*>>> -build_line_graph(const std::vector& input_polygons, SegmentLookup& segment_lookup, const std::vector& triangular_polygons) { + std::map, std::vector*>>, + std::map +> +build_line_graph(const std::vector& input_polygons, SegmentLookup& segment_lookup, const std::vector& triangular_polygons) +{ + // Build maps of triangle -> edge and edge -> triangle in order to do traversal on the 'corridor mesh' std::map, std::vector*>> segment_to_facet; std::map, std::vector*>> segment_to_input_facet; std::map, Point_2> segment_to_midpoint; std::map> midpoint_to_segment; std::map*, std::vector>> facet_to_segment; + std::map midpoint_to_edge_length; for (auto& tri : triangular_polygons) { for (size_t i = 0; i < 3; ++i) { @@ -813,6 +883,7 @@ build_line_graph(const std::vector& input_polygons, SegmentLookup& se if (p1index != input_polygons.end() && p2index != input_polygons.end() && p1index != p2index) { segment_to_midpoint[p.first] = center; midpoint_to_segment[center] = p.first; + midpoint_to_edge_length[center] = std::sqrt(CGAL::to_double(CGAL::squared_distance(p.first.first, p.first.second))); } } @@ -832,7 +903,682 @@ build_line_graph(const std::vector& input_polygons, SegmentLookup& se } } - return {line_graph, midpoint_to_segment, segment_to_input_facet}; + return {line_graph, midpoint_to_segment, segment_to_input_facet, midpoint_to_edge_length}; +} + +using DPoint = CGAL::Simple_cartesian::Point_2; +using DDir = CGAL::Simple_cartesian::Vector_2; +using DBox = std::array; + +struct CenterLineGraphData { + std::vector points; + std::vector points_double; + std::vector widths; + std::vector> edges; + std::vector> incident_edges; +}; + +struct LineRun { + Point_2 start_exact; + Point_2 end_exact; + DPoint start; + DPoint end; + DDir direction; + double avg_width; + double length; + size_t vertex_count; +}; + +struct RunBoxRecord { + size_t run_index; + DPoint start; + DPoint end; + DDir direction; + double width; + double length; + std::array corners; + DBox bbox; +}; + +struct MergedBoxRecord { + DPoint start; + DPoint end; + DDir direction; + DDir normal; + double avg_width; + double length; + size_t member_count; + std::vector members; + std::array corners; + DBox bbox; + Point_2 exact_start; + Point_2 exact_end; +}; + +struct BoxCluster { + std::vector members; + MergedBoxRecord box; +}; + +struct SnapCandidate { + size_t box_index; + double box_distance; + double line_distance; + Point_2 projection; +}; + +DDir unit(const DDir& a) { + auto n = std::sqrt(a.squared_length()); + if (n < 1.e-9) { + return {0., 0.}; + } + return a / n; +} + +DDir perpendicular(const DDir& a) { + return DDir(-a.y(), a.x()); +} + +DDir canonicalize_like(const DDir& a, const DDir& ref) { + return (a * ref) < 0. ? -a : a; +} + +DPoint to_double_point(const Point_2& p) { + return {CGAL::to_double(p.x()), CGAL::to_double(p.y())}; +} + +Point_2 to_exact_point(const DPoint& p) { + return Point_2(p.x(), p.y()); +} + +double point_line_distance(const DPoint& p, const DPoint& line_point, const DDir& line_dir) { + auto u = unit(line_dir); + auto delta = (p - line_point); + if (u.squared_length() < 1.e-18) { + return std::sqrt(delta.squared_length()); + } + return std::abs(CGAL::determinant(u.x(), u.y(), delta.x(), delta.y())); +} + +double angle_between_dirs_deg(const DDir& a, const DDir& b) { + auto u = unit(a); + auto v = unit(b); + auto c = std::abs(u * v); + if (c > 1.) { + c = 1.; + } + return std::acos(c) * 180. / 3.14159265358979323846; +} + +std::array rectangle_corners(const DPoint& start, const DPoint& end, double width) { + auto u = unit(end - start); + if (u.squared_length() < 1.e-18) { + u = {1., 0.}; + } + auto n = perpendicular(u); + auto ext = width; + auto p0 = start - u * ext; + auto p1 = end + u * ext; + auto w = n * (width / 2.); + return {p0 + w, p1 + w, p1 - w, p0 - w}; +} + +DBox aabb_from_points(const std::array& corners) { + DBox bbox{corners[0], corners[0]}; + for (auto& p : corners) { + bbox[0] = {std::min(bbox[0].x(), p.x()), std::min(bbox[0].y(), p.y())}; + bbox[1] = {std::max(bbox[1].x(), p.x()), std::max(bbox[1].y(), p.y())}; + } + return bbox; +} + +bool aabb_overlap(const DBox& a, const DBox& b, double eps = 1.e-9) { + return a[0].x() <= b[1].x() + eps && + a[1].x() + eps >= b[0].x() && + a[0].y() <= b[1].y() + eps && + a[1].y() + eps >= b[0].y(); +} + +CenterLineGraphData make_center_line_graph_data( + const std::map>& line_graph, + const std::map& midpoint_to_edge_length) +{ + CenterLineGraphData graph; + std::map point_to_index; + + auto ensure_point = [&](const Point_2& p) { + auto it = point_to_index.find(p); + if (it != point_to_index.end()) { + return it->second; + } + auto i = graph.points.size(); + point_to_index[p] = i; + graph.points.push_back(p); + graph.points_double.push_back(to_double_point(p)); + auto wt = midpoint_to_edge_length.find(p); + graph.widths.push_back(wt == midpoint_to_edge_length.end() ? 0. : wt->second); + graph.incident_edges.emplace_back(); + return i; + }; + + for (auto& p : line_graph) { + ensure_point(p.first); + for (auto& q : p.second) { + ensure_point(q); + } + } + + std::set> seen_edges; + for (auto& p : line_graph) { + auto i = ensure_point(p.first); + for (auto& q : p.second) { + auto j = ensure_point(q); + if (i == j) { + continue; + } + auto e = i < j ? std::make_pair(i, j) : std::make_pair(j, i); + if (seen_edges.insert(e).second) { + auto k = graph.edges.size(); + graph.edges.push_back(e); + graph.incident_edges[e.first].push_back(k); + graph.incident_edges[e.second].push_back(k); + } + } + } + + return graph; +} + +double segment_width(const CenterLineGraphData& graph, const std::pair& edge) { + return 0.5 * (graph.widths[edge.first] + graph.widths[edge.second]); +} + +bool edge_supports_same_line( + const DPoint& seed_a, + const DPoint& seed_b, + const DPoint& test_a, + const DPoint& test_b, + double angle_tol_deg = 3., + double line_dist_tol = 0.15) +{ + auto d_seed = seed_b - seed_a; + auto d_test = test_b - test_a; + if (d_seed.squared_length() < 1.e-18 || d_test.squared_length() < 1.e-18) { + return false; + } + if (angle_between_dirs_deg(d_seed, d_test) > angle_tol_deg) { + return false; + } + return + point_line_distance(test_a, seed_a, d_seed) <= line_dist_tol && + point_line_distance(test_b, seed_a, d_seed) <= line_dist_tol; +} + +std::vector runs_from_graph(const CenterLineGraphData& graph, double angle_tol_deg = 3., double line_dist_tol = 0.15) { + std::vector visited(graph.edges.size(), false); + std::vector runs; + + for (size_t seed_ei = 0; seed_ei < graph.edges.size(); ++seed_ei) { + if (visited[seed_ei]) { + continue; + } + + const auto& seed_edge = graph.edges[seed_ei]; + auto seed_a = graph.points_double[seed_edge.first]; + auto seed_b = graph.points_double[seed_edge.second]; + auto seed_dir = seed_b - seed_a; + if (seed_dir.squared_length() < 1.e-18) { + visited[seed_ei] = true; + continue; + } + + std::vector queue = {seed_ei}; + std::set component_edges; + + while (!queue.empty()) { + auto ei = queue.back(); + queue.pop_back(); + if (!component_edges.insert(ei).second) { + continue; + } + + const auto& edge = graph.edges[ei]; + std::array vertices = {edge.first, edge.second}; + for (auto v : vertices) { + for (auto ej : graph.incident_edges[v]) { + if (ej == ei || visited[ej] || component_edges.count(ej)) { + continue; + } + const auto& candidate = graph.edges[ej]; + auto test_a = graph.points_double[candidate.first]; + auto test_b = graph.points_double[candidate.second]; + if (edge_supports_same_line(seed_a, seed_b, test_a, test_b, angle_tol_deg, line_dist_tol)) { + queue.push_back(ej); + } + } + } + } + + for (auto ei : component_edges) { + visited[ei] = true; + } + + std::set component_vertices; + auto ref = unit(seed_dir); + DDir direction_sum{0., 0.}; + double total_length = 0.; + double weighted_width_sum = 0.; + + for (auto ei : component_edges) { + const auto& edge = graph.edges[ei]; + component_vertices.insert(edge.first); + component_vertices.insert(edge.second); + + auto d = graph.points_double[edge.second] - graph.points_double[edge.first]; + auto u = canonicalize_like(unit(d), ref); + direction_sum = direction_sum + u; + + auto len = std::sqrt(d.squared_length()); + total_length += len; + weighted_width_sum += len * segment_width(graph, edge); + } + + auto run_direction = direction_sum.squared_length() < 1.e-18 ? ref : unit(direction_sum); + + double min_t = std::numeric_limits::infinity(); + double max_t = -std::numeric_limits::infinity(); + size_t start_index = *component_vertices.begin(); + size_t end_index = start_index; + for (auto vi : component_vertices) { + auto t = (graph.points_double[vi] - CGAL::ORIGIN) * run_direction; + if (t < min_t) { + min_t = t; + start_index = vi; + } + if (t > max_t) { + max_t = t; + end_index = vi; + } + } + + auto avg_width = total_length < 1.e-9 ? segment_width(graph, seed_edge) : weighted_width_sum / total_length; + + runs.push_back({ + graph.points[start_index], + graph.points[end_index], + graph.points_double[start_index], + graph.points_double[end_index], + run_direction, + avg_width, + std::sqrt((graph.points_double[end_index] - graph.points_double[start_index]).squared_length()), + component_vertices.size() + }); + } + + return runs; +} + +std::vector build_run_box_records(const std::vector& runs) { + std::vector records; + records.reserve(runs.size()); + for (size_t i = 0; i < runs.size(); ++i) { + auto corners = rectangle_corners(runs[i].start, runs[i].end, runs[i].avg_width); + records.push_back({ + i, + runs[i].start, + runs[i].end, + unit(runs[i].end - runs[i].start), + runs[i].avg_width, + runs[i].length, + corners, + aabb_from_points(corners) + }); + } + return records; +} + +template +std::pair projected_interval_on_axis(const T& box, const DDir& axis_u) { + auto u = unit(axis_u); + auto ta = (box.start - CGAL::ORIGIN) * u; + auto tb = (box.end - CGAL::ORIGIN) * u; + return {std::min(ta, tb), std::max(ta, tb)}; +} + +double interval_overlap_length(const std::pair& a, const std::pair& b) { + return std::max(0., std::min(a.second, b.second) - std::max(a.first, b.first)); +} + +template +double boxes_overlap_along_merge_axis(const T& a, const T& b) { + auto d1 = unit(a.end - a.start); + auto d2 = unit(b.end - b.start); + if (d1 * d2 < 0.) { + d2 = {-d2.x(), -d2.y()}; + } + auto merge_axis = unit(d1 + d2); + if (merge_axis.squared_length() < 1.e-18) { + merge_axis = d1; + } + + auto i1 = projected_interval_on_axis(a, merge_axis); + auto i2 = projected_interval_on_axis(b, merge_axis); + auto overlap = interval_overlap_length(i1, i2); + auto small_length = std::min(i1.second - i1.first, i2.second - i2.first); + if (small_length < 1.e-9) { + return false; + } + return overlap / small_length; +} + +MergedBoxRecord merge_cluster_to_box(const std::vector& member_indices, const std::vector& records) { + auto ref = records[member_indices.front()].direction; + DDir direction_sum{0., 0.}; + for (auto i : member_indices) { + auto u = canonicalize_like(records[i].direction, ref); + direction_sum = direction_sum + u * std::max(records[i].length, 1.e-9); + } + + auto u = direction_sum.squared_length() < 1.e-18 ? ref : unit(direction_sum); + auto n = perpendicular(u); + + double tmin = std::numeric_limits::infinity(); + double tmax = -std::numeric_limits::infinity(); + double smin = std::numeric_limits::infinity(); + double smax = -std::numeric_limits::infinity(); + + for (auto i : member_indices) { + for (auto& corner : records[i].corners) { + auto t = (corner - CGAL::ORIGIN) * u; + auto s = (corner - CGAL::ORIGIN) * n; + tmin = std::min(tmin, t); + tmax = std::max(tmax, t); + smin = std::min(smin, s); + smax = std::max(smax, s); + } + } + + auto width = smax - smin; + auto sc = (smin + smax) / 2.; + auto start = u * tmin + n * sc; + auto end = u * tmax + n * sc; + auto corners = rectangle_corners(CGAL::ORIGIN + start, CGAL::ORIGIN + end, width); + + MergedBoxRecord box{ + CGAL::ORIGIN + start, + CGAL::ORIGIN + end, + u, + n, + width, + std::sqrt((end - start).squared_length()), + member_indices.size(), + member_indices, + corners, + aabb_from_points(corners), + to_exact_point(CGAL::ORIGIN + start), + to_exact_point(CGAL::ORIGIN + end) + }; + return box; +} + +std::pair merge_score(const MergedBoxRecord& a, const MergedBoxRecord& b) { + auto ang = angle_between_dirs_deg(a.direction, b.direction); + auto center_a = ((a.start - CGAL::ORIGIN) + (a.end - CGAL::ORIGIN)) / 2.; + auto center_b = ((b.start - CGAL::ORIGIN) + (b.end - CGAL::ORIGIN)) / 2.; + return {ang, std::sqrt((center_b - center_a).squared_length())}; +} + +bool clusters_can_merge(const BoxCluster& a, const BoxCluster& b, double angle_tol_deg = 5., double axis_overlap_ratio_limit = 0.5) { + if (!aabb_overlap(a.box.bbox, b.box.bbox)) { + return false; + } + if (angle_between_dirs_deg(a.box.direction, b.box.direction) > angle_tol_deg) { + return false; + } + if (boxes_overlap_along_merge_axis(a.box, b.box) > axis_overlap_ratio_limit) { + auto a_center = CGAL::ORIGIN + ((a.box.start - CGAL::ORIGIN) + (a.box.end - CGAL::ORIGIN)) / 2.; + auto b_center = CGAL::ORIGIN + ((b.box.start - CGAL::ORIGIN) + (b.box.end - CGAL::ORIGIN)) / 2.; + auto a_dir = a.box.direction; + auto b_dir = b.box.direction; + auto dist = a.box.length < b.box.length ? point_line_distance(a_center, b_center, b_dir) : point_line_distance(b_center, a_center, a_dir); + auto ref = a.box.length < b.box.length ? a.box.avg_width : b.box.avg_width; + return dist < (ref / 4.); + } + return true; +} + +std::vector merge_intersecting_parallel_boxes_iterative(const std::vector& runs) { + auto records = build_run_box_records(runs); + std::vector clusters; + clusters.reserve(records.size()); + for (size_t i = 0; i < records.size(); ++i) { + clusters.push_back({{i}, merge_cluster_to_box({i}, records)}); + } + + while (true) { + std::optional> best_pair; + std::pair best_score; + + for (size_t i = 0; i < clusters.size(); ++i) { + for (size_t j = i + 1; j < clusters.size(); ++j) { + if (!clusters_can_merge(clusters[i], clusters[j])) { + continue; + } + auto score = merge_score(clusters[i].box, clusters[j].box); + if (!best_pair || score < best_score) { + best_pair = std::make_pair(i, j); + best_score = score; + } + } + } + + if (!best_pair) { + break; + } + + auto i = best_pair->first; + auto j = best_pair->second; + std::vector members = clusters[i].members; + members.insert(members.end(), clusters[j].members.begin(), clusters[j].members.end()); + auto merged = BoxCluster{members, merge_cluster_to_box(members, records)}; + + std::vector next_clusters; + next_clusters.reserve(clusters.size() - 1); + for (size_t k = 0; k < clusters.size(); ++k) { + if (k != i && k != j) { + next_clusters.push_back(std::move(clusters[k])); + } + } + next_clusters.push_back(std::move(merged)); + clusters = std::move(next_clusters); + } + + std::vector merged_boxes; + merged_boxes.reserve(clusters.size()); + for (auto& cluster : clusters) { + merged_boxes.push_back(cluster.box); + } + return merged_boxes; +} + +Point_2 project_point_to_line_exact(const Point_2& p, const MergedBoxRecord& box) { + auto d = box.exact_end - box.exact_start; + if (d.squared_length() == 0) { + return box.exact_start; + } + auto t = ((p - box.exact_start) * d) / d.squared_length(); + return box.exact_start + d * t; +} + +boost::optional intersect_infinite_lines_exact(const MergedBoxRecord& a, const MergedBoxRecord& b) { + if (a.exact_start == a.exact_end || b.exact_start == b.exact_end) { + return boost::none; + } + auto x = CGAL::intersection(CGAL::Line_2(a.exact_start, a.exact_end), CGAL::Line_2(b.exact_start, b.exact_end)); + if (!x) { + return boost::none; + } + if (auto* xp = variant_get(&*x)) { + return *xp; + } + return boost::none; +} + +double point_to_oriented_box_distance(const DPoint& p, const MergedBoxRecord& box) { + auto d = box.end - box.start; + auto L = std::sqrt(d.squared_length()); + if (L < 1.e-9) { + return std::sqrt((p - box.start).squared_length()); + } + + auto u = d / L; + auto n = perpendicular(u); + auto rel = p - box.start; + auto t = rel * u; + auto s = rel * n; + + auto tmin = -box.avg_width / 2.; + auto tmax = L + box.avg_width / 2.; + auto smin = -box.avg_width / 2.; + auto smax = box.avg_width / 2.; + + double dt = 0.; + if (t < tmin) { + dt = tmin - t; + } else if (t > tmax) { + dt = t - tmax; + } + + double ds = 0.; + if (s < smin) { + ds = smin - s; + } else if (s > smax) { + ds = s - smax; + } + + return std::hypot(dt, ds); +} + +std::map> snap_points_to_box_axes( + const CenterLineGraphData& graph, + const std::vector& boxes) +{ + std::vector snapped_points(graph.points.size()); + + for (size_t i = 0; i < graph.points.size(); ++i) { + if (boxes.empty()) { + snapped_points[i] = graph.points[i]; + continue; + } + + std::vector candidates; + candidates.reserve(boxes.size()); + for (size_t j = 0; j < boxes.size(); ++j) { + candidates.push_back({ + j, + point_to_oriented_box_distance(graph.points_double[i], boxes[j]), + point_line_distance(graph.points_double[i], boxes[j].start, boxes[j].direction), + project_point_to_line_exact(graph.points[i], boxes[j]) + }); + } + + std::vector containing; + for (auto& candidate : candidates) { + if (candidate.box_distance <= 1.e-9) { + containing.push_back(candidate); + } + } + + auto less = [](const SnapCandidate& a, const SnapCandidate& b) { + if (a.line_distance != b.line_distance) { + return a.line_distance < b.line_distance; + } + return a.box_distance < b.box_distance; + }; + + if (containing.size() >= 2) { + std::sort(containing.begin(), containing.end(), less); + auto& c1 = containing[0]; + auto& c2 = containing[1]; + if (angle_between_dirs_deg(boxes[c1.box_index].direction, boxes[c2.box_index].direction) > 8.) { + if (auto x = intersect_infinite_lines_exact(boxes[c1.box_index], boxes[c2.box_index])) { + snapped_points[i] = *x; + continue; + } + } + snapped_points[i] = c1.projection; + continue; + } + + if (containing.size() == 1) { + snapped_points[i] = containing[0].projection; + continue; + } + + auto best = *std::min_element(candidates.begin(), candidates.end(), [](const SnapCandidate& a, const SnapCandidate& b) { + if (a.box_distance != b.box_distance) { + return a.box_distance < b.box_distance; + } + return a.line_distance < b.line_distance; + }); + snapped_points[i] = best.projection; + } + + std::map> adjacency; + for (auto& edge : graph.edges) { + auto a = snapped_points[edge.first]; + auto b = snapped_points[edge.second]; + if (a == b) { + continue; + } + adjacency[a].insert(b); + adjacency[b].insert(a); + } + + std::map> snapped_graph; + for (auto& p : adjacency) { + snapped_graph[p.first] = {p.second.begin(), p.second.end()}; + } + return snapped_graph; +} + +Graph2D join_segment_runs( + DebugWriter& debug, + const std::map>& line_graph, + const std::map& midpoint_to_edge_length) +{ + auto graph = make_center_line_graph_data(line_graph, midpoint_to_edge_length); + auto runs = runs_from_graph(graph); + runs.erase(std::remove_if(runs.begin(), runs.end(), [](const LineRun& run) { + return run.vertex_count <= 5; + }), runs.end()); + + std::vector run_polygons; + for (auto& r : runs) { + auto ps = rectangle_corners(r.start, r.end, r.avg_width); + std::array exact_corners; + std::transform(ps.begin(), ps.end(), exact_corners.begin(), [](const DPoint& p) { + return to_exact_point(p); + }); + run_polygons.emplace_back(exact_corners.begin(), exact_corners.end()); + } + debug.write_polygons(run_polygons, "initial_runs"); + run_polygons.clear(); + + auto boxes = merge_intersecting_parallel_boxes_iterative(runs); + + for (auto& r : boxes) { + auto ps = rectangle_corners(r.start, r.end, r.avg_width); + std::array exact_corners; + std::transform(ps.begin(), ps.end(), exact_corners.begin(), [](const DPoint& p) { + return to_exact_point(p); + }); + run_polygons.emplace_back(exact_corners.begin(), exact_corners.end()); + } + debug.write_polygons(run_polygons, "merged_boxes"); + + auto snapped_graph = snap_points_to_box_axes(graph, boxes); + return Graph2D(snapped_graph); } std::set> find_triangles(const std::map>& line_graph) { @@ -1118,66 +1864,88 @@ std::list> extend_end_vertices_based_on_input( const Graph2D& G, const std::map>& midpoint_to_segment, const std::map, std::vector*>>& segment_to_input_facet, - const Polygon_list& inner_offset, - const SegmentLookup& segment_lookup + const Polygon_list& outer_perimiter, + const SegmentLookup& segment_lookup, + const K::FT& max_projection_distance ){ std::list> constructed_segments; - for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) { - if (it->second.size() == 1) { - auto& M = it->first; + std::set processed_vertices; - const std::pair* q = nullptr; + while (true) { + // The idea was to peal off 1-degree vertices when projecting them did not result into + // nearby intersections with the outer perimiter. This in case there would be turns near + // the perimeter, which would be eliminated by pealing off the vertices, which would then + // require out of the loop because of invalidated iterators. For now we decided to stick + // to a projection of the vertex onto the perimeter segment when the projection distance + // exceeds a threshold. + bool broke_out = false; - if (midpoint_to_segment.find(M) == midpoint_to_segment.end()) { - typename K::FT min_sq_distance = std::numeric_limits::infinity(); - for (auto& pa : midpoint_to_segment) { - if (CGAL::squared_distance(pa.first, M) < min_sq_distance) { - q = &pa.second; - min_sq_distance = CGAL::squared_distance(pa.first, M); - } + for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) { + if (it->second.size() == 1) { + auto& M = it->first; + + if (processed_vertices.find(M) != processed_vertices.end()) { + continue; } - } else { - q = &midpoint_to_segment.find(M)->second; - } - if (q == nullptr) { - continue; - } + const std::pair* q = nullptr; - bool handled_as_graph_path = false; + if (midpoint_to_segment.find(M) == midpoint_to_segment.end()) { + typename K::FT min_sq_distance = std::numeric_limits::infinity(); + for (auto& pa : midpoint_to_segment) { + if (CGAL::squared_distance(pa.first, M) < min_sq_distance) { + q = &pa.second; + min_sq_distance = CGAL::squared_distance(pa.first, M); + } + } + } else { + q = &midpoint_to_segment.find(M)->second; + } - // distance from unioned - shoot ray? - 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 - CGAL::Ray_2 ray(incoming, M - incoming); - // intersect ray with boundary - boost::optional> closest_segment; - boost::optional> closest_intersection_point; - K::FT sq_distance_along_ray = std::numeric_limits::infinity(); - for (auto jt = bnd.edges_begin(); jt != bnd.edges_end(); ++jt) { - const auto& seg = *jt; - auto x = CGAL::intersection(ray, seg); - if (x) { - if (auto* xp = variant_get>(&*x)) { - auto dist = ((*xp) - M).squared_length(); - if (dist < sq_distance_along_ray) { - closest_segment = seg; - closest_intersection_point = *xp; - sq_distance_along_ray = dist; + if (q == nullptr) { + continue; + } + + bool handled_as_graph_path = false; + + // distance from unioned - shoot ray? + if (segment_to_input_facet.find(*q)->second.size() == 2) { + for (auto& bnd : outer_perimiter) { + // 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 + CGAL::Ray_2 ray(incoming, M - incoming); + + // intersect ray with boundary + boost::optional> closest_segment; + boost::optional> closest_intersection_point; + K::FT sq_distance_along_ray = std::numeric_limits::infinity(); + for (auto jt = bnd.edges_begin(); jt != bnd.edges_end(); ++jt) { + const auto& seg = *jt; + auto x = CGAL::intersection(ray, seg); + if (x) { + if (auto* xp = variant_get>(&*x)) { + auto dist = ((*xp) - M).squared_length(); + if (dist < sq_distance_along_ray) { + if (dist < (max_projection_distance * max_projection_distance)) { + closest_segment = seg; + closest_intersection_point = *xp; + sq_distance_along_ray = dist; + } else { + + } + } } } } - } - if (closest_intersection_point) { - constructed_segments.push_front({M, *closest_intersection_point}); - break; + if (closest_intersection_point) { + constructed_segments.push_front({M, *closest_intersection_point}); + processed_vertices.insert(M); + break; #if 0 Graph2D GGG(bnd); GGG.refine(*GGG.query(*closest_intersection_point, 0.01), *closest_intersection_point); @@ -1217,12 +1985,35 @@ std::list> extend_end_vertices_based_on_input( break; } #endif - } else { - std::cerr << "Warning: no intersection found when extending end vertex, this will likely result in invalid topology" << std::endl; + } else { + + // Loop over boundary segments, and project point onto it, take the closest + K::FT closest_distance = std::numeric_limits::infinity(); + boost::optional> closest_point; + for (auto& poly : outer_perimiter) { + for (auto jt = poly.edges_begin(); jt != poly.edges_end(); ++jt) { + auto seg = *jt; + auto Pp = seg.supporting_line().projection(M); + if (seg.has_on(Pp)) { + auto d = CGAL::squared_distance(Pp, M); + if (d < (max_projection_distance * max_projection_distance)) { + if (d < closest_distance) { + closest_distance = d; + closest_point = Pp; + } + } + } + } + } + + if (closest_point) { + constructed_segments.push_front({M, *closest_point}); + processed_vertices.insert(M); + } + } } } } - } #if 0 if (!handled_as_graph_path) { @@ -1267,6 +2058,11 @@ std::list> extend_end_vertices_based_on_input( constructed_segments.push_front({avg, R}); } #endif + } + } + + if (!broke_out) { + break; } } @@ -1355,8 +2151,6 @@ 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 { @@ -1367,7 +2161,112 @@ class Segment_2_less { } }; -void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { +std::vector arrangement_cell_iou(Arrangement_2& left, Arrangement_2& right) { + + using Walk_pl = CGAL::Arr_walk_along_line_point_location; + Walk_pl walk_pl(right); + + std::set visited_faces_on_right; + + std::vector return_values; + + for (auto it = left.faces_begin(); it != left.faces_end(); ++it) { + if (!it->is_unbounded()) { + // convert arr facet to polygon with holes + auto polygon_exterior = circ_to_poly(it->outer_ccb()); + Polygon_with_holes_2 pwh(polygon_exterior); + for (auto hit = it->inner_ccbs_begin(); hit != it->inner_ccbs_end(); ++hit) { + pwh.add_hole(circ_to_poly(*hit)); + } + + CGAL::Polygon_triangulation_decomposition_2 decompositor; + std::vector temp; + decompositor(pwh, std::back_inserter(temp)); + + std::set visited_points; + + while (true) { + // select triangle edge that has largest squared edge length times distance from polygon exterior + K::FT max_score = -std::numeric_limits::infinity(); + Point_2 best_point; + for (auto& tri : temp) { + for (size_t i = 0; i < 3; ++i) { + size_t j = (i + 1) % 3; + auto& pi = tri.vertex(i); + auto& pj = tri.vertex(j); + + auto center_point = CGAL::ORIGIN + (((pi - CGAL::ORIGIN) + (pj - CGAL::ORIGIN)) / 2); + + K::FT min_dist = std::numeric_limits::infinity(); + for (auto eit = polygon_exterior.edges_begin(); eit != polygon_exterior.edges_end(); ++eit) { + auto ep = eit->source(); + auto eq = eit->target(); + Segment_2 seg(ep, eq); + auto dist = CGAL::squared_distance(center_point, seg); + if (dist < min_dist) { + min_dist = dist; + } + } + + auto sq_length = CGAL::squared_distance(pi, pj); + + auto score = sq_length * min_dist; + if (score > max_score && visited_points.count(center_point) == 0) { + max_score = score; + best_point = center_point; + } + } + } + + auto res = walk_pl.locate(best_point); + if (auto* v = variant_get(&res)) { + if (visited_faces_on_right.count(*v) > 0) { + return_values.push_back(0); + } else { + // convert arr facet to polygon with holes + auto polygon_exterior = circ_to_poly((*v)->outer_ccb()); + Polygon_with_holes_2 pwh_right(polygon_exterior); + for (auto hit = (*v)->inner_ccbs_begin(); hit != (*v)->inner_ccbs_end(); ++hit) { + pwh_right.add_hole(circ_to_poly(*hit)); + } + + // compute intersection over union of pwh and the original polygon + if (CGAL::do_intersect(pwh, pwh_right)) { + std::vector result; + CGAL::intersection(pwh, pwh_right, std::back_inserter(result)); + typename K::FT intersection_area = 0; + for (auto& r : result) { + auto poly_area = r.outer_boundary().area(); + for (auto& h : r.holes()) { + poly_area -= h.area(); + } + intersection_area += poly_area; + } + CGAL::Polygon_with_holes_2 poly12; + CGAL::join(pwh, pwh_right, poly12); + typename K::FT union_area = poly12.outer_boundary().area(); + for (auto& h : poly12.holes()) { + union_area -= h.area(); + } + return_values.push_back(intersection_area / union_area); + } else { + return_values.push_back(0); + } + } + visited_faces_on_right.insert(*v); + break; + } else { + // Not in facet on right, retry another point + continue; + } + } + } + } + + return return_values; +} + +void clean_noisy_paths(DebugWriter& debug_output, Arrangement_2& arr, SegmentLookup& segment_lookup, double& threshold) { using SK = CGAL::Simple_cartesian; CGAL::Cartesian_converter C{}; @@ -1418,7 +2317,7 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { auto [dv, dl] = get_dir(s); best = std::min(best, angle(dv)); } - return (best + 0.1) / own_length; + return (best + 0.01) / own_length; }; std::map badnesses; @@ -1426,7 +2325,6 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { badnesses[e] = edge_badness(e); } - double thr; { std::vector tmp; tmp.reserve(badnesses.size()); @@ -1435,12 +2333,12 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { } std::nth_element(tmp.begin(), tmp.begin() + tmp.size() / 2, tmp.end()); double med = tmp[tmp.size() / 2]; - thr = 10.0 * med; + threshold = 4.0 * med; } std::set bad_edges; for (auto& p : badnesses) { - if (p.second > thr) { + if (p.second > threshold) { bad_edges.insert(p.first); } } @@ -1568,7 +2466,46 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { return best_x; }; + auto process_modifications = [&]( + Arrangement_2& arr_, + const std::set>& to_remove_, + const std::vector>& to_insert_) { + 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)); + } + }; + + size_t path_index = 0; for (auto& path : bad_paths) { + decltype(to_remove) to_remove_this_path; + decltype(to_insert) to_insert_this_path; + + for (size_t i = 0; i < path.size() - 1; ++i) { + auto& a = path[i]; + auto& b = path[i + 1]; + + debug_output.write_segment(a, b, "arr_bad_path path_nr_" + std::to_string(path_index)); + } + auto x = collapse_path(path); if (!x) { // std::cerr << "Unable to collapse path, skipping" << std::endl; @@ -1589,12 +2526,10 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { 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; + // std::cerr << "Collapsing path would create near-duplicate vert to previous path, skipping" << std::endl; continue; } } @@ -1604,75 +2539,314 @@ void clean_noisy_paths(Arrangement_2& arr, SegmentLookup& segment_lookup) { auto& b = path[i + 1]; if (a < b) { to_remove.insert({a, b}); + to_remove_this_path.insert({a, b}); } else { to_remove.insert({b, a}); + to_remove_this_path.insert({b, a}); } } auto s = path.front(); auto t = path.back(); if (s != *x) { to_insert.push_back({s, *x}); + to_insert_this_path.push_back({s, *x}); + + debug_output.write_segment(s, *x, "corrected_path path_nr_" + std::to_string(path_index)); } if (t != *x) { to_insert.push_back({t, *x}); + to_insert_this_path.push_back({t, *x}); + + debug_output.write_segment(t, *x, "corrected_path path_nr_" + std::to_string(path_index)); } + + path_index += 1; + +#if 1 + process_modifications(arr, to_remove_this_path, to_insert_this_path); +#else + auto arr_copy = arr; + process_modifications(arr_copy, to_remove_this_path, to_insert_this_path); + auto ious = arrangement_cell_iou(arr, arr_copy); + for (auto& iou : ious) { + std::cerr << " - cell iou: " << CGAL::to_double(iou) << std::endl; + } + std::swap(arr_copy, arr); +#endif } - /* - using Walk_pl = CGAL::Arr_walk_along_line_point_location; - Walk_pl walk_pl(arr); + process_modifications(arr, to_remove, to_insert); +} - 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); +template +void next_circular(typename Vec::const_iterator& it, const Vec& vec) { + std::advance(it, 1); + if (it == vec.end()) { + it = vec.begin(); + } +} +template +void previous_circular(typename Vec::const_iterator& it, const Vec& vec) { + if (it == vec.begin()) { + it = vec.end(); + } + std::advance(it, -1); +} - if ((*v)->point() != e.first) { - std::cerr << "Warning: unable to locate vertex for edge removal, skipping" << std::endl; - continue; +template +std::size_t circular_distance(typename Vec::const_iterator first, + typename Vec::const_iterator last, + const Vec& vec) { + if (first <= last) { + return static_cast(last - first); + } + return static_cast(vec.end() - first) + static_cast(last - vec.begin()); +} + +template +std::pair +longest_wrapping_true_run(const Vec& v, Pred pred) { + using It = typename Vec::const_iterator; + + const auto n = v.size(); + if (n == 0) { + return {v.end(), v.end()}; + } + + // Find best non-wrapping run + std::size_t best_len = 0; + std::size_t best_start = 0; + + std::size_t curr_len = 0; + std::size_t curr_start = 0; + + for (std::size_t i = 0; i < n; ++i) { + if (pred(v[i])) { + if (curr_len == 0) { + curr_start = i; + } + ++curr_len; + if (curr_len > best_len) { + best_len = curr_len; + best_start = curr_start; } } else { - std::cerr << "Warning: unable to locate vertex for edge removal, skipping" << std::endl; - continue; + curr_len = 0; } } - */ - 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; + // Count leading true + std::size_t leading = 0; + while (leading < n && pred(v[leading])) { + ++leading; + } + + // All true + if (leading == n) { + return {v.begin(), v.end()}; + } + + // Count trailing true + std::size_t trailing = 0; + while (trailing < n && pred(v[n - 1 - trailing])) { + ++trailing; + } + + // Wrapped run = [n - trailing, n) + [0, leading) + const std::size_t wrapped_len = leading + trailing; + + if (wrapped_len > best_len) { + It first = v.begin() + static_cast(n - trailing); + It last = v.begin() + static_cast(leading); + return {first, last}; + } + + It first = v.begin() + static_cast(best_start); + It last = first + static_cast(best_len); + return {first, last}; +} + +void clean_noisy_bounds(DebugWriter& debug_output, Arrangement_2& arr, SegmentLookup& segment_lookup, double threshold) { + 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(); + }; + + 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.01) / own_length; + }; + + size_t facet_index = 0; + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it, ++facet_index) { + if (!it->is_unbounded()) { + std::set> to_remove; + std::vector> to_insert; + + std::vector segs; + std::vector vertices; + std::vector halfedges; + + auto circ = it->outer_ccb(); + do { + auto a = circ->source()->point(); + auto b = circ->target()->point(); + segs.emplace_back(a, b); + vertices.push_back(circ->source()); + halfedges.push_back(circ); + ++circ; + } while (circ != it->outer_ccb()); + + std::vector badnesses; + for (auto& e : segs) { + badnesses.push_back(edge_badness(e)); + } + + auto bit = std::min_element(badnesses.begin(), badnesses.end()); + if (*bit > threshold) { + // std::cerr << "All edges are good, skipping" << std::endl; + continue; + } + + auto it_pair = longest_wrapping_true_run(badnesses, [&](double d) { return d > threshold; }); + auto N = circular_distance(it_pair.first, it_pair.second, badnesses); + + if (N == 0) { + // std::cerr << "Unable to find run of bad edges, skipping" << std::endl; + continue; + } + + std::vector> incoming_paths; + + auto jt = it_pair.first; + for (std::size_t k = 0; k < N; ++k, next_circular(jt, badnesses)) { + + auto he = halfedges[std::distance(badnesses.cbegin(), jt)]; + to_remove.insert({he->source()->point(), he->target()->point()}); + debug_output.write_segment(he->source()->point(), he->target()->point(), "arr_bad_bound facet_" + std::to_string(facet_index)); + + Arrangement_2::Vertex_handle v = he->source(); + + // circle around other edges onto v + Arrangement_2::Halfedge_around_vertex_circulator first, curr; + first = curr = v->incident_halfedges(); + do { + Arrangement_2::Vertex_handle u = curr->source(); + if (curr->face() != it && curr->twin()->face() != it) { + + // loop until we find a 3-degree vertex, or we come back to the start + std::vector path{v->point(), u->point()}; + auto he = curr; + + while (u->degree() == 2 && u != v && path.size() < 10) { + std::vector hes; + + { + Arrangement_2::Halfedge_around_vertex_circulator first, curr; + first = curr = u->incident_halfedges(); + do { + hes.push_back(curr); + curr++; + } while (curr != first); + } + + auto next_he = hes.front() != he && hes.front() != he->twin() ? hes.front() : hes.back(); + auto next_v = next_he->target() != u ? next_he->target() : next_he->source(); + + path.push_back(next_v->point()); + u = next_v; + } + incoming_paths.push_back(std::move(path)); + } + } while (++curr != first); + } + + const std::size_t start = + static_cast(std::distance(badnesses.cbegin(), it_pair.first)); + + auto n = badnesses.size(); + + auto wrap = [n](std::ptrdiff_t i) -> std::size_t { + i %= static_cast(n); + if (i < 0) { + i += static_cast(n); + } + return static_cast(i); + }; + + const std::size_t ib = start; + const std::size_t ia = wrap(static_cast(start) - 1); + const std::size_t ic = wrap(static_cast(start + N)); + const std::size_t id = wrap(static_cast(start + N + 1)); + + auto a = vertices.begin() + static_cast(ia); + auto b = vertices.begin() + static_cast(ib); + auto c = vertices.begin() + static_cast(ic); + auto d = vertices.begin() + static_cast(id); + + CGAL::Ray_2 r1((*a)->point(), (*b)->point()); + CGAL::Ray_2 r2((*d)->point(), (*c)->point()); + + auto x = CGAL::intersection(r1, r2); + if (x) { + if (auto* xp = variant_get>(&*x)) { + to_insert.emplace_back((*b)->point(), *xp); + to_insert.emplace_back((*c)->point(), *xp); + + debug_output.write_segment((*b)->point(), *xp, "corrected_bound facet_" + std::to_string(facet_index)); + debug_output.write_segment((*c)->point(), *xp, "corrected_bound facet_" + std::to_string(facet_index)); + } + } else { + CGAL::Line_2 r1((*a)->point(), (*b)->point()); + CGAL::Line_2 r2((*d)->point(), (*c)->point()); + + auto x = CGAL::intersection(r1, r2); + if (x) { + if (auto* xp = variant_get>(&*x)) { + to_insert.emplace_back((*b)->point(), *xp); + to_insert.emplace_back((*c)->point(), *xp); + + debug_output.write_segment((*b)->point(), *xp, "corrected_bound facet_" + std::to_string(facet_index)); + debug_output.write_segment((*c)->point(), *xp, "corrected_bound facet_" + std::to_string(facet_index)); + } + } } } - 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) { @@ -1725,20 +2899,31 @@ class timer { public: class entry { public: + entry() {} + entry(std::map::const_iterator start_it) : start_it(start_it) {} + void stop() { - auto end = std::chrono::high_resolution_clock::now(); - auto duration = std::chrono::duration(end - start_it->second).count(); - std::cerr << "Timing for " << start_it->first << ": " << duration << " ms" << std::endl; + if (start_it) { + auto end = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration(end - start_it.value()->second).count(); + std::cerr << "Timing for " << start_it.value()->first << ": " << duration << " ms" << std::endl; + } } private: - std::map::const_iterator start_it; + std::optional::const_iterator> start_it; }; + timer(bool enabled = true) : enabled_(enabled) {} + entry start(const std::string& name) { - return entry(timings_.insert({name, std::chrono::high_resolution_clock::now()}).first); + if (enabled_) { + return entry(timings_.insert({name, std::chrono::high_resolution_clock::now()}).first); + } else { + return entry(); + } } private: @@ -1746,27 +2931,30 @@ class timer { std::string, std::chrono::high_resolution_clock::time_point> timings_; + + bool enabled_; }; -void arrange_cgal_polygons(const std::vector& input_polygons_, std::vector& output_polygons, double polygon_offset_distance = -1.) { +void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std::vector& input_polygons_, std::vector& output_polygons, double polygon_offset_distance = -1.) { 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; -#ifdef SVGFILL_DEBUG - auto t = std::time(nullptr); - auto tm = *std::localtime(&t); + DebugWriter debug_output; + if (settings.debug_output) { + auto t = std::time(nullptr); + auto tm = *std::localtime(&t); - std::ostringstream oss; - oss << std::put_time(&tm, "arrangement_%Y%m%d%H%M%S"); - auto now = oss.str(); - DebugWriter debug_output(true, now); -#else - DebugWriter debug_output(false, ""); -#endif + std::ostringstream oss; + oss << std::put_time(&tm, "arrangement_%Y%m%d%H%M%S"); + auto now = oss.str(); + debug_output = DebugWriter(true, now); + } else { + debug_output = DebugWriter(false, ""); + } - timer timer; + timer timer(settings.debug_output); auto t0 = timer.start("input"); @@ -1794,7 +2982,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v t0.stop(); t0 = timer.start("overlap elimination"); - eliminate_overlaps(OVERLAP_RESOLUTION_DISTANCE, input_polygons); + eliminate_overlaps(debug_output, OVERLAP_RESOLUTION_DISTANCE, input_polygons); t0.stop(); @@ -1813,79 +3001,79 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v debug_output.write_polygons(input_polygons, "processed_input"); -#if 1 - t0 = timer.start("outer perimeter"); - - // Find the outer perimeter using offset - union - negative offset - std::vector offset_polygons; - for (auto& r : input_polygons) { - auto R = r; - if (!R.is_counterclockwise_oriented()) { - R.reverse_orientation(); - } - - // Overlap removal can also result in close points causing problems when converted into non-exact nt - remove_close_points(R); - - auto ps = create_and_convert_offset_polygon(polygon_offset_distance, R); - for (auto& p : ps) { - if (!p.is_simple()) { - throw std::runtime_error("Complex polygon originated from offset"); - } - } - offset_polygons.insert(offset_polygons.end(), ps.begin(), ps.end()); - } - - debug_output.write_polygons(offset_polygons, "offset_input"); - - // Perform Boolean union on the offset polygons - std::vector unioned_polygons; - CGAL::join(offset_polygons.begin(), offset_polygons.end(), std::back_inserter(unioned_polygons)); - - if (unioned_polygons.size() > 1) { - // @todo this is currently one of the major limitations in the code that still can be eliminated - // by grouping the input polygons by their perimiter polygon in unioned_polygons - std::sort(unioned_polygons.begin(), unioned_polygons.end(), [](auto& p, auto& q) { return p.outer_boundary().area() > q.outer_boundary().area(); }); - } - - debug_output.write_polygon(unioned_polygons.front().outer_boundary(), "offset_joined"); - - Polygon_2 fused_removed_close_points = unioned_polygons.front().outer_boundary(); - remove_close_points(fused_removed_close_points, 1.e-4); - - // Apply negative offset to get the outer perimeter 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()), - - // 1.e-8 even was too little and still resulted in slivers of triangle around the perimeter - -polygon_offset_distance - OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT, - fused_removed_close_points); - - 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 + if (settings.outer_perimiter_algo == 0) { + t0 = timer.start("outer perimeter"); + + // Find the outer perimeter using offset - union - negative offset + std::vector offset_polygons; + for (auto& r : input_polygons) { + auto R = r; + if (!R.is_counterclockwise_oriented()) { + R.reverse_orientation(); + } + + // Overlap removal can also result in close points causing problems when converted into non-exact nt + remove_close_points(R); + + auto ps = create_and_convert_offset_polygon(polygon_offset_distance, R); + for (auto& p : ps) { + if (!p.is_simple()) { + throw std::runtime_error("Complex polygon originated from offset"); + } + } + offset_polygons.insert(offset_polygons.end(), ps.begin(), ps.end()); + } + + debug_output.write_polygons(offset_polygons, "offset_input"); + + // Perform Boolean union on the offset polygons + std::vector unioned_polygons; + CGAL::join(offset_polygons.begin(), offset_polygons.end(), std::back_inserter(unioned_polygons)); + + if (unioned_polygons.size() > 1) { + // @todo this is currently one of the major limitations in the code that still can be eliminated + // by grouping the input polygons by their perimiter polygon in unioned_polygons + std::sort(unioned_polygons.begin(), unioned_polygons.end(), [](auto& p, auto& q) { return p.outer_boundary().area() > q.outer_boundary().area(); }); + } + + debug_output.write_polygon(unioned_polygons.front().outer_boundary(), "offset_joined"); + + Polygon_2 fused_removed_close_points = unioned_polygons.front().outer_boundary(); + remove_close_points(fused_removed_close_points, 1.e-4); + + // Apply negative offset to get the outer perimeter polygon + 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()), + + // 1.e-8 even was too little and still resulted in slivers of triangle around the perimeter + -polygon_offset_distance - OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT, + fused_removed_close_points); + + 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); + } + outer_perimiter.emplace_back(cycle.begin(), cycle.end()); + } t0.stop(); t0 = timer.start("corridor creation"); @@ -1911,8 +3099,10 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // subdivide difference_result to have better more detailed triangulation and therefore less-pronounced artefacts in midpoint network + auto subdivision_length = polygon_offset_distance / settings.subdivision_factor; + for (auto& pwh : difference_result) { - difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 8., pwh)); + difference_result_subdivided.push_back(subdivide_polygon(subdivision_length, pwh)); // difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 64., pwh)); } @@ -1940,7 +3130,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v SegmentLookup segment_lookup(input_polygons); - auto [line_graph, midpoint_to_segment, segment_to_input_facet] = build_line_graph(input_polygons, segment_lookup, triangular_polygons); + auto [line_graph, midpoint_to_segment, segment_to_input_facet, midpoint_to_edge_length] = build_line_graph(input_polygons, segment_lookup, triangular_polygons); for (auto& p : line_graph) { for (auto& q : p.second) { debug_output.write_segment(p.first, q, "network_1"); @@ -1950,38 +3140,48 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v t0.stop(); t0 = timer.start("center line cleaning"); - - auto triangles = find_triangles(line_graph); - // For every triangle found in the network we eliminate one edge to break the cycle - // The edge we eliminate is the edge with the greatest angle with any of it's neighbours - auto eliminated_segments = eliminate_triangles(line_graph); + Graph2D G; + if (settings.line_cleaning_algo == 0) { + G = join_segment_runs(debug_output, line_graph, midpoint_to_edge_length); + Arrangement_2 arr; + G.to_arrangement(arr); + Graph2D G2; + G2.from_arrangement(arr); + eliminate_colinear_vertices(G2); + G = G2; + for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { + debug_output.write_segment(it->first, it->second, "network_2"); + } + } else { + auto eliminated_segments = eliminate_triangles(line_graph); - Graph2D G2(line_graph); - for (auto& e : eliminated_segments) { - debug_output.write_segment(e.first, e.second, "eliminated"); - G2.remove_edge(e.first, e.second); - } + Graph2D G2(line_graph); + for (auto& e : eliminated_segments) { + debug_output.write_segment(e.first, e.second, "eliminated"); + G2.remove_edge(e.first, e.second); + } - auto G = G2.weld_vertices(); + G = G2.weld_vertices(); - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - debug_output.write_segment(it->first, it->second, "network_2"); - } + for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { + debug_output.write_segment(it->first, it->second, "network_2"); + } - eliminate_colinear_vertices(G); + eliminate_colinear_vertices(G); - edge_slide(G); + edge_slide(G); - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - debug_output.write_segment(it->first, it->second, "network_3"); + for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { + debug_output.write_segment(it->first, it->second, "network_3"); + } } t0.stop(); t0 = timer.start("topology"); - auto segments = extend_end_vertices_based_on_input(G, midpoint_to_segment, segment_to_input_facet, outer_perimiter, segment_lookup); + auto segments = extend_end_vertices_based_on_input(G, midpoint_to_segment, segment_to_input_facet, outer_perimiter, segment_lookup, subdivision_length * 4); // Now plot the edges on an arrangement in order to find planar cycles // and merge the corridor-halves with their neighbouring input polygon @@ -1997,31 +3197,31 @@ 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(); - if (poly.vertex(i) == poly.vertex(j)) { - continue; + if (settings.topology_reconstruction_algo != 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(); + if (poly.vertex(i) == poly.vertex(j)) { + continue; + } + 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)); } - 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; @@ -2047,12 +3247,17 @@ 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 + if (settings.topology_reconstruction_algo != 0) { + fuse_corridor_halves_with_input(arr, G, segment_lookup, input_polygons, debug_output); + } + + if (settings.perform_cleanup && settings.line_cleaning_algo != 0) { + remove_colinear_vertices(arr); + double threshold; + clean_noisy_paths(debug_output, arr, segment_lookup, threshold); + remove_colinear_vertices(arr); + clean_noisy_bounds(debug_output, arr, segment_lookup, threshold); + } t0.stop(); @@ -2068,8 +3273,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v #ifndef SVGFILL_MAIN -bool svgfill::arrange_polygons(const std::vector& polygons, std::vector& arranged) -{ +bool svgfill::arrange_polygons(arrange_polygon_settings settings, const std::vector& polygons, std::vector& arranged) { std::vector cgal_polygons, cgal_polygons_out; std::transform(polygons.begin(), polygons.end(), std::back_inserter(cgal_polygons), [](auto& poly) { Polygon_2 result; @@ -2078,7 +3282,7 @@ bool svgfill::arrange_polygons(const std::vector& polygons, }); return result; }); - arrange_cgal_polygons(cgal_polygons, cgal_polygons_out); + arrange_cgal_polygons(settings, cgal_polygons, cgal_polygons_out); std::transform(cgal_polygons_out.begin(), cgal_polygons_out.end(), std::back_inserter(arranged), [](auto& poly) { svgfill::polygon_2 result; std::transform(poly.begin(), poly.end(), std::back_inserter(result.boundary), [](auto& pt) { @@ -2128,7 +3332,7 @@ int main(int argc, char** argv) { input_polygons.back().push_back(CGAL::Point_2(x, y)); } } - arrange_cgal_polygons(input_polygons, output); + arrange_cgal_polygons(arrange_polygon_settings{}, input_polygons, output); break; } return 0; @@ -2141,7 +3345,7 @@ int main(int argc, char** argv) { input_polygons = { rect1, rect2, rect3, rect4, rect5 }; } - arrange_cgal_polygons(input_polygons, output); + arrange_cgal_polygons(arrange_polygon_settings{}, input_polygons, output); return 0; } diff --git a/src/svgfill/src/graph_2d.h b/src/svgfill/src/graph_2d.h index da2b5014ec..d19418ff62 100644 --- a/src/svgfill/src/graph_2d.h +++ b/src/svgfill/src/graph_2d.h @@ -346,6 +346,13 @@ public: } } + template + void from_arrangement(T& arr) { + for (auto it = arr.edges_begin(); it != arr.edges_end(); ++it) { + insert(it->source()->point(), it->target()->point()); + } + } + void assert_symmetric() { #ifdef SVGFILL_DEBUG #if 0 diff --git a/src/svgfill/src/svgfill.cpp b/src/svgfill/src/svgfill.cpp index 8a2a1bb008..cf45881c93 100644 --- a/src/svgfill/src/svgfill.cpp +++ b/src/svgfill/src/svgfill.cpp @@ -483,6 +483,18 @@ public: return ps; } + size_t delete_same_facet_edge_pairs() { + size_t n_deleted = 0; + for (auto it = arr.edges_begin(); it != arr.edges_end();) { + decltype(it) current = it++; + if (current->face() == current->twin()->face()) { + arr.remove_edge(current); + n_deleted++; + } + } + return n_deleted; + } + void merge(const std::vector& edge_indices) { if (edge_indices.empty()) { return; diff --git a/src/svgfill/src/svgfill.h b/src/svgfill/src/svgfill.h index 396fc9c924..2588636a8d 100644 --- a/src/svgfill/src/svgfill.h +++ b/src/svgfill/src/svgfill.h @@ -67,6 +67,7 @@ namespace svgfill { virtual std::vector get_face_pairs() = 0; virtual size_t num_edges() = 0; virtual size_t num_faces() = 0; + virtual size_t delete_same_facet_edge_pairs() = 0; }; class SVGFILL_API context { @@ -101,6 +102,7 @@ namespace svgfill { void write(std::vector>&); size_t num_edges() { return arr_->num_edges(); } size_t num_faces() { return arr_->num_faces(); } + size_t delete_same_facet_edge_pairs() { return arr_->delete_same_facet_edge_pairs(); } ~context() { delete arr_; @@ -113,7 +115,25 @@ namespace svgfill { SVGFILL_API std::string polygons_to_svg(const std::vector>& polygons, bool random_color=false); SVGFILL_API std::string polygons_to_svg(const std::vector& polygons, bool random_color = false); SVGFILL_API bool svg_to_polygons(const std::string& data, const boost::optional& class_name, std::vector& polygons); - SVGFILL_API bool arrange_polygons(const std::vector& polygons, std::vector& arranged); -} + + struct SVGFILL_API arrange_polygon_settings { + bool debug_output = false; + // -1: compute from average edge length + double polygon_offset_distance = -1.; + // 0: use offset - union - negative offset to find the outer perimeter + // 1: radial walk along vertices; exact, but can only reuse vertices, not create new positions by means of intersections + int outer_perimiter_algo = 0; + // 0: outer perimiter and corridor center lines + // 1: input polygons, corridor center lines and segments connecting corridor center lines to input polygons + int topology_reconstruction_algo = 0; + // 0: join segment runs + // 1: local badness reduction + int line_cleaning_algo = 0; + bool perform_cleanup = true; + double subdivision_factor = 16.; + }; + + SVGFILL_API bool arrange_polygons(arrange_polygon_settings settings, const std::vector& polygons, std::vector& arranged); + } #endif