diff --git a/src/svgfill/src/arrange_polygons.cpp b/src/svgfill/src/arrange_polygons.cpp index 7bddc2e41a..0271d76146 100644 --- a/src/svgfill/src/arrange_polygons.cpp +++ b/src/svgfill/src/arrange_polygons.cpp @@ -164,31 +164,6 @@ boost::optional subtract_retain_largest(const T& lhs, const T& rhs) { return boost::none; } -// Function to write polygons as line segments in OBJ format -void write_polygon_to_obj(std::ofstream& ofs, size_t& vertex_index, bool as_line, const Polygon_2& polygon, const std::string& name) { - ofs << "o " << name << "\n"; // Object name - - // Write vertices - for (auto vit = polygon.vertices_begin(); vit != polygon.vertices_end(); ++vit) { - ofs << "v " << CGAL::to_double(vit->x()) << " " << CGAL::to_double(vit->y()) << " 0\n"; - } - - if (as_line) { - // Write line segments (edges) - for (size_t j = 0; j < polygon.size(); ++j) { - ofs << "l " << vertex_index + j << " " << vertex_index + (j + 1) % polygon.size() << "\n"; - } - } else { - ofs << "f"; - for (size_t j = 0; j < polygon.size(); ++j) { - ofs << " " << vertex_index + j; - } - ofs << "\n"; - } - - vertex_index += polygon.size(); -} - Polygon_2 circ_to_poly(typename Arrangement_2::Ccb_halfedge_const_circulator circ) { Polygon_2 poly; @@ -208,27 +183,6 @@ Polygon_with_holes_2 circ_to_poly(typename Arrangement_2::Ccb_halfedge_const_cir return poly; } -void write_polygon_to_svg(std::ostream& ofs, const Polygon_2& polygon) { - ofs << "x()) << "," << CGAL::to_double(vit->y()) << " "; - } - ofs << "\" style=\"fill:none;stroke-width:1\" />\n"; -} - -// Function to write a Polygon_with_holes_2 to an SVG file -void write_polygon_with_holes_to_svg(std::ostream& ofs, const Polygon_with_holes_2& polygon_with_holes) { - // Write the outer boundary (main polygon) - if (!polygon_with_holes.is_unbounded()) { - write_polygon_to_svg(ofs, polygon_with_holes.outer_boundary()); - } - - // Write the holes (if any) with a different color (e.g., red) - for (auto hit = polygon_with_holes.holes_begin(); hit != polygon_with_holes.holes_end(); ++hit) { - write_polygon_to_svg(ofs, *hit); - } -} - Polygon_2 fuse_with_offset(const std::vector& polygons, double polygon_offset_distance) { // Find the outer perimeter using offset - union - negative offset std::vector offset_polygons; @@ -285,78 +239,77 @@ Polygon_2 fuse_with_offset(const std::vector& polygons, double polygo return inner_offset.front(); } -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; - // 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; +double estimate_polygon_offset_distance(const std::vector& polygons) { + double total_edge_length = 0.; + size_t num_edges = 0; + for (auto& p : polygons) { + for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { + total_edge_length += std::sqrt(CGAL::to_double(CGAL::squared_distance(it->start(), it->end()))); + num_edges += 1; + } + } + return total_edge_length / num_edges / 2; +} - if (polygon_offset_distance < 0.) { - double total_edge_length = 0.; - size_t num_edges = 0; - for (auto& p : input_polygons_) { - for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { - total_edge_length += std::sqrt(CGAL::to_double(CGAL::squared_distance(it->start(), it->end()))); - num_edges += 1; +void clean_polygon(Polygon_2& poly) { + // Ensure counterclockwise orientation and remove duplicate last point if present also remove close points + if (!poly.is_counterclockwise_oriented()) { + poly.reverse_orientation(); + } + std::vector> ps(poly.begin(), poly.end()); + if (ps.front() == ps.back()) { + ps.pop_back(); + } + poly = Polygon_2(ps.begin(), ps.end()); + remove_close_points(poly); +} + +void smooth_polygon(double factor, Polygon_2& poly) { + auto ps = create_and_convert_offset_polygon(-factor, poly); + if (ps.size() == 1) { + auto r2 = ps.front(); + ps = create_and_convert_offset_polygon(+factor, r2); + if (ps.size() == 1) { + poly = ps.front(); + } + } +} + +template +void split_self_intersecting_polygon(const CGAL::Polygon_2& poly, OutIt output_it) { + if (poly.is_simple()) { + *output_it++ = poly; + return; + } + Arrangement_2 arr; + for (auto it = poly.edges_begin(); it != poly.edges_end(); ++it) { + CGAL::insert(arr, Segment_2(it->start(), it->end())); + } + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (it->is_unbounded()) { + for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { + auto inner = circ_to_poly(*jt); + // reverse because it's an inner bound to the infinite outer facet + inner.reverse_orientation(); + *output_it++ = inner; } } - polygon_offset_distance = total_edge_length / num_edges / 2; } +} - auto input_polygons__ = input_polygons_; - decltype(input_polygons__) input_polygons; - - for (auto& i : input_polygons__) { - std::vector> ps(i.begin(), i.end()); - if (ps.front() == ps.back()) { - ps.pop_back(); - } - input_polygons.emplace_back(ps.begin(), ps.end()); - } - - for (auto& polygon : input_polygons) { - if (!polygon.is_counterclockwise_oriented()) { - polygon.reverse_orientation(); - } - } - - for (auto& polygon : input_polygons) { - remove_close_points(polygon); - } - -#ifdef SVGFILL_DEBUG - std::ofstream obj("obj.obj"); - size_t vi = 1; - - std::ofstream svg("svg.svg"); - svg << "\n"; - - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "input_poly_" + std::to_string(std::distance(input_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } - - obj << std::flush; -#endif - +std::set> +find_overlaps(const std::vector& polygons) { typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + std::vector boxes; + std::vector>> input_triangulated; - std::vector boxes; - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons.begin(); it != polygons.end(); ++it) { constexpr double offset = 1.e-3; auto b = it->bbox(); boxes.emplace_back( CGAL::Bbox_2(b.xmin() - offset, b.ymin() - offset, b.xmax() + offset, b.ymax() + offset), - std::distance(input_polygons.begin(), it) - ); - - if (!it->is_simple()) { -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, *it, "self-intersecting"); -#endif - throw std::runtime_error("Self-intersecting input"); - } + std::distance(polygons.begin(), it)); CGAL::Polygon_triangulation_decomposition_2 decompositor; std::vector temp; @@ -378,10 +331,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v bool registered_overlap = false; for (auto& t2 : input_triangulated[b.handle()]) { if (CGAL::squared_distance(t1, t2) < (1.e-3 * 1.e-3)) { - overlaps.insert({ - (a.handle() < b.handle()) ? a.handle() : b.handle(), - (a.handle() < b.handle()) ? b.handle() : a.handle() - }); + overlaps.insert({(a.handle() < b.handle()) ? a.handle() : b.handle(), + (a.handle() < b.handle()) ? b.handle() : a.handle()}); registered_overlap = true; break; } @@ -393,297 +344,319 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } }); - if (true) { - // solve overlaps by means of subtraction - // loop over overlaps and subtract the smaller polygon from the larger one + return overlaps; +} - std::set eliminated_polies; - std::map overlap_counts; - for (auto& p : overlaps) { - overlap_counts[p.first]++; - overlap_counts[p.second]++; +class DebugWriter { + public: + DebugWriter(bool enabled, const std::string& filename_prefix) + : enabled_(enabled) { + if (enabled_) { + obj.open(filename_prefix + ".obj"); + vi = 1; + svg.open(filename_prefix + ".svg"); + svg << "\n"; } - - for (const auto& edge : overlaps) { - // Skip eliminated - if (eliminated_polies.find(edge.first) != eliminated_polies.end() || - eliminated_polies.find(edge.second) != eliminated_polies.end()) { - continue; - } - - // Many overlaps indicate an aggregated polygon, skip them - /* - if (overlap_counts[edge.first] > 10 || overlap_counts[edge.second] > 10) { - if (overlap_counts[edge.first] > 10) { - eliminated_polies.insert(edge.first); - } - if (overlap_counts[edge.second] > 10) { - eliminated_polies.insert(edge.second); - } - continue; - } - */ - - // these are pointers now, because otherwise swap would not work? - auto* poly1 = &input_polygons[edge.first]; - auto* poly2 = &input_polygons[edge.second]; - - // Populate eliminated_polies with small polygons - // This can happen over time when modifications are made to the polygons to solve overlaps - bool skip = false; - if (poly1->area() < 1.e-2) { - eliminated_polies.insert(edge.first); - skip = true; - } - if (poly2->area() < 1.e-2) { - eliminated_polies.insert(edge.second); - skip = true; - } - // Small slivers are also just eliminated - if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly1))) { - eliminated_polies.insert(edge.first); - skip = true; - } - if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly2))) { - eliminated_polies.insert(edge.second); - skip = true; - } - if (skip) { - continue; - } - - // Skip polygons that have a very high intersection over union - // ratio, which indicates that they are very likely duplicates - if (CGAL::do_intersect(*poly1, *poly2)) { - std::vector result; - CGAL::intersection(*poly1, *poly2, 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(*poly1, *poly2, poly12); - typename K::FT union_area = poly12.outer_boundary().area(); - for (auto& h : poly12.holes()) { - union_area -= h.area(); - } - if (union_area > 0 && intersection_area / union_area > 0.99) { - // std::cerr << intersection_area / union_area << std::endl; - eliminated_polies.insert(edge.first); - continue; - } - } - - if (!(poly1->is_simple() && poly2->is_simple())) { - continue; - } - - { - std::vector result; - // std::cerr << poly1.area() << " " << poly2.area() << std::endl; - // std::cerr.flush(); - - boost::optional mp1, mp2, mp3, mp4; - bool swap = false; - - swap = poly1->area() <= poly2->area(); - if (swap) { - std::swap(poly1, poly2); - } - - bool success = false; - if ((mp1 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE, *poly2)))) { - if ((mp2 = subtract_retain_largest(*poly1, *mp1))) { - if ((mp3 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE * 2, *mp2)))) { - if ((mp4 = subtract_retain_largest(*poly2, *mp3))) { - *poly1 = *mp2; - *poly2 = *mp4; - success = true; - } - } - } - } - - /* - if (swap) { - // swap back to retain original ordering - // what's the point in swapping back here? - std::swap(poly1, poly2); - } - */ - - if (!success) { - eliminated_polies.insert(swap ? edge.first : edge.second); - continue; - } - } - } - - // iterate over the eliminated polygons and remove them from the input polygons - for (auto it = eliminated_polies.rbegin(); it != eliminated_polies.rend(); ++it) { - input_polygons.erase(input_polygons.begin() + *it); + } + ~DebugWriter() { + if (enabled_) { + svg << "\n"; + obj << std::flush; + obj.close(); + svg.close(); } } + void write_polygon(const Polygon_2& polygon, const std::string& name) { + if (enabled_) { + write_polygon_to_obj_(obj, vi, true, polygon, name); + write_polygon_to_svg_(svg, polygon, name); + obj << std::flush; + } + } + + void write_segment(const Point_2& p, const Point_2& q, const std::string& name) { + if (enabled_) { + if (last_segment_name_ != name) { + last_segment_name_ = name; + obj << "o " << name << "\n"; + } + obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; + obj << "v " << CGAL::to_double(q.x()) << " " << CGAL::to_double(q.y()) << " 0\n"; + obj << "l " << vi++; + obj << " " << vi++ << "\n"; + + svg << ""; + + obj << std::flush; + } + } + + void write_polygon(const Polygon_with_holes_2& polygon, const std::string& name) { + if (enabled_) { + write_polygon(polygon.outer_boundary(), name); + for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) { + write_polygon(*hit, name); + } + } + } + + void write_polygons(const std::vector& polygons, const std::string& name) { + if (enabled_) { + size_t i = 0; + for (auto& polygon : polygons) { + write_polygon_to_obj_(obj, vi, true, polygon, name + "_" + std::to_string(i++)); + write_polygon_to_svg_(svg, polygon, name); + } + obj << std::flush; + } + } + + void write_polygons(const std::vector& polygons, const std::string& name) { + if (enabled_) { + size_t i = 0; + for (auto& polygon : polygons) { + write_polygon_to_obj_(obj, vi, true, polygon.outer_boundary(), name + "_" + std::to_string(i)); + write_polygon_to_svg_(svg, polygon.outer_boundary(), name); + for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) { + write_polygon_to_obj_(obj, vi, true, *hit, name + "_" + std::to_string(i)); + write_polygon_to_svg_(svg, *hit, name); + } + } + obj << std::flush; + } + } + + private: + std::ofstream obj; + size_t vi; + std::ofstream svg; + bool enabled_; + std::string last_segment_name_; + + void write_polygon_to_svg_(std::ostream& ofs, const Polygon_2& polygon, const std::string& class_name = "") { + ofs << "x()) << "," << -CGAL::to_double(vit->y()) << " "; + } + ofs << "\"/>\n"; + } + + void write_polygon_to_obj_(std::ofstream& ofs, size_t& vertex_index, bool as_line, const Polygon_2& polygon, const std::string& name) { + ofs << "o " << name << "\n"; // Object name + + // Write vertices + for (auto vit = polygon.vertices_begin(); vit != polygon.vertices_end(); ++vit) { + ofs << "v " << CGAL::to_double(vit->x()) << " " << CGAL::to_double(vit->y()) << " 0\n"; + } + + if (as_line) { + // Write line segments (edges) + for (size_t j = 0; j < polygon.size(); ++j) { + ofs << "l " << vertex_index + j << " " << vertex_index + (j + 1) % polygon.size() << "\n"; + } + } else { + ofs << "f"; + for (size_t j = 0; j < polygon.size(); ++j) { + ofs << " " << vertex_index + j; + } + ofs << "\n"; + } + + vertex_index += polygon.size(); + } +}; + +void eliminate_overlaps(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 + + std::set eliminated_polies; + /* - if constexpr (false) { - // solve overlap by means of union into components - std::vector> adj(input_polygons.size()); - for (const auto& edge : overlaps) { - adj[edge.first].push_back(edge.second); - adj[edge.second].push_back(edge.first); - } - - std::vector visited(input_polygons.size(), false); - std::vector> connected_components; - - for (size_t v = 0; v < input_polygons.size(); ++v) { - if (!visited[v]) { - connected_components.emplace_back(); - - std::stack stack; - stack.push(v); - visited[v] = true; - - while (!stack.empty()) { - size_t u = stack.top(); - stack.pop(); - connected_components.back().push_back(u); - - for (size_t neighbor : adj[u]) { - if (!visited[neighbor]) { - visited[neighbor] = true; - stack.push(neighbor); - } - } - } - } - } - - std::vector fused_polies; - - for (auto& comp : connected_components) { - std::vector comp_polies; - if (comp.size() == 1) { - fused_polies.push_back(input_polygons[comp.front()]); - } else { - for (auto& c : comp) { - comp_polies.push_back(input_polygons[c]); - } - fused_polies.push_back(fuse_with_offset(comp_polies, 1.e-2)); - } - } - -#ifdef SVGFILL_DEBUG - for (auto it = fused_polies.begin(); it != fused_polies.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "fused_poly_" + std::to_string(std::distance(fused_polies.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - - input_polygons = fused_polies; + std::map overlap_counts; + for (auto& p : overlaps) { + overlap_counts[p.first]++; + overlap_counts[p.second]++; } */ - { - // [NB Nov 6] we cannot do this anymore because it could revert the spacing between input polygons - // that touch in the corner. - // Now that overlaps/touches at corners are handled more locally only a small indent is produced - // which would be undone by means of an inset+offset. - // - // [NB Nov 10] this is actually still necessary though, but we apply a much smaller distance now - // to keep the overlap eliminations in tact - // - // Inset-offset to remove tiny details that may cause enourmous spikes in offsets - for (auto& r : input_polygons) { - auto ps = create_and_convert_offset_polygon(-polygon_offset_distance / 10000., r); - if (ps.size() == 1) { - auto r2 = ps.front(); - ps = create_and_convert_offset_polygon(+polygon_offset_distance / 10000., r2); - if (ps.size() == 1) { - r = ps.front(); + auto overlaps = find_overlaps(polygons); + + for (const auto& edge : overlaps) { + // Skip eliminated + if (eliminated_polies.find(edge.first) != eliminated_polies.end() || + eliminated_polies.find(edge.second) != eliminated_polies.end()) { + continue; + } + + // Many overlaps indicate an aggregated polygon, skip them + /* + if (overlap_counts[edge.first] > 10 || overlap_counts[edge.second] > 10) { + if (overlap_counts[edge.first] > 10) { + eliminated_polies.insert(edge.first); + } + if (overlap_counts[edge.second] > 10) { + eliminated_polies.insert(edge.second); + } + continue; + } + */ + + // these are pointers now, because otherwise swap would not work? + auto* poly1 = &polygons[edge.first]; + auto* poly2 = &polygons[edge.second]; + + // @todo this is applied during overlap processing, maybe better after the boolean operation, + // because they can be come small or narrow when overlaps are resolved + + // Populate eliminated_polies with small polygons + // This can happen over time when modifications are made to the polygons to solve overlaps + bool skip = false; + if (poly1->area() < 1.e-2) { + eliminated_polies.insert(edge.first); + skip = true; + } + if (poly2->area() < 1.e-2) { + eliminated_polies.insert(edge.second); + skip = true; + } + // Small slivers are also just eliminated + if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly1))) { + eliminated_polies.insert(edge.first); + skip = true; + } + if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly2))) { + eliminated_polies.insert(edge.second); + skip = true; + } + if (skip) { + continue; + } + + // Skip polygons that have a very high intersection over union + // ratio, which indicates that they are very likely duplicates + if (CGAL::do_intersect(*poly1, *poly2)) { + std::vector result; + CGAL::intersection(*poly1, *poly2, 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(*poly1, *poly2, poly12); + typename K::FT union_area = poly12.outer_boundary().area(); + for (auto& h : poly12.holes()) { + union_area -= h.area(); + } + if (union_area > 0 && intersection_area / union_area > 0.99) { + // std::cerr << intersection_area / union_area << std::endl; + eliminated_polies.insert(edge.first); + continue; + } + } + + if (!(poly1->is_simple() && poly2->is_simple())) { + continue; + } + + { + std::vector result; + + boost::optional mp1, mp2, mp3, mp4; + bool swap = false; + + swap = poly1->area() <= poly2->area(); + if (swap) { + std::swap(poly1, poly2); + } + + bool success = false; + if ((mp1 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE, *poly2)))) { + if ((mp2 = subtract_retain_largest(*poly1, *mp1))) { + if ((mp3 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE * 2, *mp2)))) { + if ((mp4 = subtract_retain_largest(*poly2, *mp3))) { + *poly1 = *mp2; + *poly2 = *mp4; + success = true; + } + } + } + } + + if (!success) { + eliminated_polies.insert(swap ? edge.first : edge.second); + continue; } } } -#ifdef SVGFILL_DEBUG - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "processed_input_poly_" + std::to_string(std::distance(input_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); + // iterate over the eliminated polygons and remove them from the input polygons + for (auto it = eliminated_polies.rbegin(); it != eliminated_polies.rend(); ++it) { + polygons.erase(polygons.begin() + *it); } -#endif +} - // Unfortunately CGAL does not seem to have a ready to use aabb primitive for segments in 2D, - // so we have to use 3D segments and aabb tree for 2D polygons. - std::list> all_segs; - std::unordered_map*, decltype(input_polygons.begin())> seg_to_poly; +class SegmentLookup { + public: + typedef std::vector::const_iterator PolygonIt; - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - for (auto eit = it->edges_begin(); eit != it->edges_end(); ++eit) { - CGAL::Segment_3 seg3d( - CGAL::Point_3(eit->source().x(), eit->source().y(), 0), - CGAL::Point_3(eit->target().x(), eit->target().y(), 0) - ); - all_segs.push_back(seg3d); - seg_to_poly[&all_segs.back()] = it; - } - } - - using TreeTraits = CGAL::AABB_traits>::iterator>>; - using Tree = CGAL::AABB_tree; - - Tree tree(all_segs.begin(), all_segs.end()); - tree.accelerate_distance_queries(); - - auto input_polygon_boundary = - [&](const Point_2& p, double tol = 1e-5) -> decltype(input_polygons.begin()) + SegmentLookup(const std::vector& polygons) + : polygons_ref_(polygons) { + // Unfortunately CGAL does not seem to have a ready to use aabb primitive for segments in 2D, + // so we have to use 3D segments and aabb tree for 2D polygons. + for (auto it = polygons.begin(); it != polygons.end(); ++it) { + for (auto eit = it->edges_begin(); eit != it->edges_end(); ++eit) { + CGAL::Segment_3 seg3d( + CGAL::Point_3(eit->source().x(), eit->source().y(), 0), + CGAL::Point_3(eit->target().x(), eit->target().y(), 0)); + all_segs.push_back(seg3d); + seg_to_poly[&all_segs.back()] = it; + } + } + tree_ = Tree(all_segs.begin(), all_segs.end()); + tree_.accelerate_distance_queries(); + } + + // This part is the most computationally expensive. Caching effectively halves the lookup time here, since every vertex on the subdivided corridor mesh has on average two outgoing edges. + PolygonIt input_polygon_boundary(const Point_2& p, double tol = 1e-5) { + auto it = input_polygon_boundary_cache_.find(p); + if (it != input_polygon_boundary_cache_.end()) { + return it->second; + } + // Find closest point & corresponding segment - auto closest = tree.closest_point_and_primitive(CGAL::Point_3(p.x(), p.y(), 0)); + auto closest = tree_.closest_point_and_primitive(CGAL::Point_3(p.x(), p.y(), 0)); const auto& closest_pt = closest.first; auto seg_ptr = &*closest.second; double d = CGAL::to_double(CGAL::squared_distance(p, Point_2(closest_pt.x(), closest_pt.y()))); + + PolygonIt res; if (d < (tol * tol)) { - return seg_to_poly.find(seg_ptr)->second; + res = seg_to_poly.find(seg_ptr)->second; + } else { + res = polygons_ref_.end(); } - return input_polygons.end(); + + input_polygon_boundary_cache_[p] = res; + return res; }; - /* - auto input_polygon_boundary = [&input_polygons](const CGAL::Point_2& p, double tol = 1.e-5) { - // unfortunately some imprecision slept into the code so we can't - // so we can't just use has_on_boundary() anymore - double D = std::numeric_limits::infinity(); - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - for (auto jt = it->edges_begin(); jt != it->edges_end(); ++jt) { - const auto& seg = *jt; - auto d = std::sqrt(CGAL::to_double(CGAL::squared_distance(seg, p))); - if (d < D) { - D = d; - } - if (d < tol) { - return it; - } - } - } - return input_polygons.end(); - }; - */ - - auto close_input_point = [&input_polygons](const CGAL::Point_2& P) { + std::pair> close_input_point(const CGAL::Point_2& P) const { + // @todo use tree CGAL::Point_2 closest; double closest_distance = std::numeric_limits::infinity(); - auto input_it = input_polygons.end(); + auto input_it = polygons_ref_.end(); // unfortunately some imprecision slept into the code so we can't // so we can't just use has_on_boundary() anymore - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons_ref_.begin(); it != polygons_ref_.end(); ++it) { for (auto& p : *it) { auto d = std::sqrt(CGAL::to_double(CGAL::squared_distance(P, p))); if (d < closest_distance) { @@ -697,14 +670,16 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v return std::make_pair(input_it, closest); }; - auto project_input_point = [&input_polygons](const CGAL::Point_2& P) { + std::pair> project_input_point(const CGAL::Point_2& P) const { + // @todo use tree + CGAL::Point_2 closest; typename K::FT closest_sq_distance = std::numeric_limits::infinity(); - auto input_it = input_polygons.end(); + auto input_it = polygons_ref_.end(); // unfortunately some imprecision slept into the code so we can't // so we can't just use has_on_boundary() anymore - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons_ref_.begin(); it != polygons_ref_.end(); ++it) { for (auto jt = it->edges_begin(); jt != it->edges_end(); ++jt) { auto Pp = jt->supporting_line().projection(P); auto d = CGAL::squared_distance(Pp, P); @@ -719,224 +694,52 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v return std::make_pair(input_it, closest); }; - // 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(); - } +private: + using TreeTraits = CGAL::AABB_traits>::iterator>>; + using Tree = CGAL::AABB_tree; - // Overlap removal can also result in close points causing problems when converted into non-exact nt - remove_close_points(R); + const std::vector& polygons_ref_; + std::list> all_segs; + std::unordered_map*, PolygonIt> seg_to_poly; + Tree tree_; - auto ps = create_and_convert_offset_polygon(polygon_offset_distance, R); - for (auto& p : ps) { - if (!p.is_simple()) { - /*{ - std::cerr << "input ["; - bool first = true; - for (auto& pp : r) { - if (!first) { - std::cerr << ","; - } - first = false; - std::cerr << "(" << pp.x() << "," << pp.y() << ")"; - } - std::cerr << "]" << std::endl; - } + std::map::const_iterator> input_polygon_boundary_cache_; +}; - { - std::cerr << "["; - bool first = true; - for (auto& pp : p) { - if (!first) { - std::cerr << ","; - } - first = false; - std::cerr << "(" << pp.x() << "," << pp.y() << ")"; - } - std::cerr << "]" << std::endl; - }*/ - - throw std::runtime_error("Complex polygon originated from offset"); - } - } - offset_polygons.insert(offset_polygons.end(), ps.begin(), ps.end()); - } - -#ifdef SVGFILL_DEBUG - for (auto it = offset_polygons.begin(); it != offset_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "offset_poly_" + std::to_string(std::distance(offset_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - - // 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(); }); - } - -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, unioned_polygons.front().outer_boundary(), "offset_poly_joined"); - write_polygon_to_svg(svg, unioned_polygons.front().outer_boundary()); - -#endif - - Polygon_2 fused_removed_close_points; - { - std::vector> ps; - auto& p = unioned_polygons.front().outer_boundary(); - ps.reserve(p.size()); - auto I = p.begin(); - auto J = I + 1; - for (;; ++J) { - bool last = false; - if (J == p.end()) { - J = p.begin(); - last = true; - } - // if (CGAL::squared_distance(*I, *J) > (polygon_offset_distance * polygon_offset_distance)) { - if (CGAL::squared_distance(*I, *J) > (1.e-4 * 1.e-4)) { - ps.push_back(*J); - I = J; - } - if (last) { - break; - } - } - fused_removed_close_points = Polygon_2(ps.begin(), ps.end()); - } - - // Apply negative offset to get the outer perimeter polygon - auto inner_offset = 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); - -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, inner_offset.front(), "joined_inset"); - write_polygon_to_svg(svg, inner_offset.front()); -#endif - - /* - // there is non-insignificant chance that around the outer boundary, vertices are located in - // between of the input polyhedra, but intermediate vertices result in triangles that will no longer - // span between the two spaces with two edges and therefore cause the topological centre line - // to no run up to the center. Eliminate all vertices that are not on the polyhedral boundary of polygon. - - // this theory proved to be false. once we have topological end points in our graph that are - // connected to input polyhedra to form closed cells, we move those topological end points to - // the average of the input polyhedra corner points, thus effectively also moving them outwards. - { - for (auto& i : inner_offset) { - std::vector> ps; - for (auto& p : i) { - if (input_polygon_boundary(p, 1.e-3) != input_polygons.end()) { - ps.push_back(p); - } - } - i = Polygon_2(ps.begin(), ps.end()); +Polygon_2 subdivide_polygon(double max_distance, const Polygon_2 & p) { + std::vector points; + for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { + const auto& seg = *it; + auto num_splits = (int)std::ceil(std::sqrt(CGAL::to_double(seg.squared_length())) / max_distance) - 1; + points.push_back(seg.source()); + for (auto i = 0; i < num_splits; ++i) { + auto d = (seg.target() - seg.source()) / (num_splits + 1) * (i + 1); + points.push_back(seg.source() + d); } } + return Polygon_2(points.begin(), points.end()); +}; -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, inner_offset.front(), "joined_inset_cleaned"); - write_polygon_to_svg(svg, inner_offset.front()); -#endif - */ - - // Subtract original polygons from outer perimeter - std::vector difference_result, difference_result_subdivided; - for (auto& i : inner_offset) { - std::vector working_copy; - working_copy.emplace_back(i); - - for (auto& r : input_polygons) { - std::vector temp_working_copy; - for (auto& wc : working_copy) { - CGAL::difference(wc, r, std::back_inserter(temp_working_copy)); - } - working_copy = temp_working_copy; - } - difference_result.insert(difference_result.end(), working_copy.begin(), working_copy.end()); +Polygon_with_holes_2 subdivide_polygon(double max_distance, const Polygon_with_holes_2& pwh) { + Polygon_2 outer = subdivide_polygon(max_distance, pwh.outer_boundary()); + std::vector holes; + for (auto hit = pwh.holes_begin(); hit != pwh.holes_end(); ++hit) { + holes.push_back(subdivide_polygon(max_distance, *hit)); } + return Polygon_with_holes_2(outer, holes.begin(), holes.end()); +}; - // subdivide difference_result to have better behave triangulation - - { - const double max_distance = polygon_offset_distance / 8.; - auto subdivide_polygon = [max_distance](const Polygon_2& p) { - std::vector points; - for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { - const auto& seg = *it; - auto num_splits = (int)std::ceil(std::sqrt(CGAL::to_double(seg.squared_length())) / max_distance) - 1; - points.push_back(seg.source()); - for (auto i = 0; i < num_splits; ++i) { - auto d = (seg.target() - seg.source()) / (num_splits + 1) * (i + 1); - points.push_back(seg.source() + d); - } - } - return Polygon_2(points.begin(), points.end()); - }; - - for (auto& pwh : difference_result) { - // Subdivide outer boundary - Polygon_2 outer = subdivide_polygon(pwh.outer_boundary()); - // Subdivide holes - std::vector holes; - for (auto hit = pwh.holes_begin(); hit != pwh.holes_end(); ++hit) { - holes.push_back(subdivide_polygon(*hit)); - } - // Construct new Polygon_with_holes_2 - difference_result_subdivided.push_back(Polygon_with_holes_2(outer, holes.begin(), holes.end())); - } - } - -#ifdef SVGFILL_DEBUG - for (auto it = difference_result_subdivided.begin(); it != difference_result_subdivided.end(); ++it) { - auto i = std::distance(difference_result_subdivided.begin(), it); - write_polygon_to_obj(obj, vi, true, it->outer_boundary(), "difference_result_subdivided_" + std::to_string(i)); - write_polygon_to_svg(svg, it->outer_boundary()); - for (auto& p : it->holes()) { - write_polygon_to_obj(obj, vi, true, p, "difference_result_subdivided_" + std::to_string(i)); - write_polygon_to_svg(svg, p); - } - } -#endif - - std::list> triangular_polygons; - - for (auto& pwh : difference_result_subdivided) { - CGAL::Polygon_triangulation_decomposition_2 decompositor; - decompositor(pwh, std::back_inserter(triangular_polygons)); - } - - triangular_polygons.erase(std::remove_if(triangular_polygons.begin(), triangular_polygons.end(), [](const CGAL::Polygon_2& p) { - return CGAL::to_double(p.area()) < 1.e-8; - }), triangular_polygons.end()); - -#ifdef SVGFILL_DEBUG - for (auto it = triangular_polygons.begin(); it != triangular_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, false, *it, "tri_" + std::to_string(std::distance(triangular_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - +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) { // 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, 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*, std::vector>> facet_to_segment; for (auto& tri : triangular_polygons) { for (size_t i = 0; i < 3; ++i) { @@ -950,26 +753,14 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } - // This part is the most computationally expensive. Caching effectively halves the lookup time here, since every vertex has two outgoing edges. - std::map input_polygon_boundary_cache; - auto cached_input_polygon_boundary = [&](const Point_2& p, double tol = 1e-5) -> decltype(input_polygons.begin()) - { - auto it = input_polygon_boundary_cache.find(p); - if (it == input_polygon_boundary_cache.end()) { - auto index = input_polygon_boundary(p, tol); - input_polygon_boundary_cache[p] = index; - return index; - } else { - return it->second; - } - }; + // @todo The smarter thing to do probably after creating the corridor mesh, register segments wrt to originating input polygon(s) and maintain that mapping when subdividing // Register midpoints on the edges within the 'corridor mesh' that span multiple input polygons for (auto& p : segment_to_facet) { auto center = CGAL::ORIGIN + (((p.first.first - CGAL::ORIGIN) + (p.first.second - CGAL::ORIGIN)) / 2); - auto p1index = cached_input_polygon_boundary(p.first.first); - auto p2index = cached_input_polygon_boundary(p.first.second); + auto p1index = segment_lookup.input_polygon_boundary(p.first.first); + auto p2index = segment_lookup.input_polygon_boundary(p.first.second); segment_to_input_facet[p.first].push_back(&*p1index); segment_to_input_facet[p.first].push_back(&*p2index); @@ -978,17 +769,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v segment_to_midpoint[p.first] = center; midpoint_to_segment[center] = p.first; } - - if (p1index != input_polygons.end() && p2index != input_polygons.end() && p1index != p2index) { - segment_to_midpoint[p.first] = center; - midpoint_to_segment[center] = p.first; - } } -#ifdef SVGFILL_DEBUG - obj << "o network_1\n"; -#endif - // Observe corridor mesh topology to join edge midpoints into a network std::map> line_graph; for (auto& p : segment_to_midpoint) { @@ -1000,20 +782,15 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v decltype(segment_to_midpoint)::const_iterator it; if ((it = segment_to_midpoint.find(r)) != segment_to_midpoint.end()) { line_graph[p.second].push_back(it->second); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(p.second.x()) << " " << CGAL::to_double(p.second.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - - svg << "second.x()) << "\" y2=\"" << CGAL::to_double(it->second.y()) << "\" />"; -#endif } } } } + return {line_graph, midpoint_to_segment, segment_to_input_facet}; +} + +std::set> find_triangles(const std::map>& line_graph) { // Find triangles in this network often occuring at junctions in the corridor mesh std::set> triangles; std::function&)> find_triangles_recursive; @@ -1024,7 +801,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v const std::vector& neighbors_current = line_graph.at(path.back()); if (std::find(neighbors_current.begin(), neighbors_current.end(), path.front()) != neighbors_current.end()) { // We found a triangle, add it to the set - Triangle triangle = { path[0], path[1], path[2] }; + Triangle triangle = {path[0], path[1], path[2]}; std::sort(triangle.begin(), triangle.end()); triangles.insert(triangle); } @@ -1037,29 +814,28 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (std::find(path.begin(), path.end(), neighbor) == path.end()) { path.push_back(neighbor); find_triangles_recursive(path); - path.pop_back(); // Backtrack + path.pop_back(); // Backtrack } } }; for (auto& p : line_graph) { - std::vector ps = { p.first }; + std::vector ps = {p.first}; find_triangles_recursive(ps); } - // 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 + return triangles; +} - // non exact time, we need sqrt +std::set> eliminate_triangles(const std::map>& line_graph) { + auto triangles = find_triangles(line_graph); + + // @todo this currently uses a simple cartesian kernel for performance for support of sqrt, but + // this should be possible to rewrite as ratios/slopes in the exact kernel as well using SK = CGAL::Simple_cartesian; CGAL::Cartesian_converter C{}; -#ifdef SVGFILL_DEBUG - obj << "o eliminated\n"; -#endif - std::set> eliminated_segments; - for (auto& t : triangles) { Triangle st; std::transform(t.begin(), t.end(), st.begin(), C); @@ -1075,7 +851,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v double max_abs_dot = 0.; { - auto& ni = line_graph[t[i]]; + auto& ni = line_graph.find(t[i])->second; for (auto& n : ni) { if (std::find(t.begin(), t.end(), n) == t.end()) { // not contained in triangle @@ -1092,7 +868,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } { - auto& nj = line_graph[t[j]]; + auto& nj = line_graph.find(t[j])->second; for (auto& n : nj) { if (std::find(t.begin(), t.end(), n) == t.end()) { // not contained in triangle @@ -1106,7 +882,6 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } - } if (max_abs_dot < global_min_abs_dot) { @@ -1119,176 +894,132 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto i = global_min_abs_dot_index; auto j = (i + 2) % 3; - eliminated_segments.insert({ t[i], t[j] }); - eliminated_segments.insert({ t[j], t[i] }); - -#ifdef SVGFILL_DEBUG - obj << "v " << st[j].x() << " " << st[j].y() << " 0\n"; - obj << "v " << st[i].x() << " " << st[i].y() << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - - svg << ""; -#endif + eliminated_segments.insert({t[i], t[j]}); + eliminated_segments.insert({t[j], t[i]}); } - } - Graph2D G2(line_graph); - for (auto& e : eliminated_segments) { - G2.remove_edge(e.first, e.second); + return eliminated_segments; +} + +bool is_parallel_2degree_node(Graph2D::vertex_const_iterator vit) { + auto it = vit->second.begin(); + auto& P = *it++; + auto& Q = *it++; + auto e1 = P - vit->first; + auto e2 = vit->first - Q; + if (e1.squared_length() == 0 || e2.squared_length() == 0) { + // @todo why does this happen? + return false; } + e1 /= std::sqrt(CGAL::to_double(e1.squared_length())); + e2 /= std::sqrt(CGAL::to_double(e2.squared_length())); + return std::abs(CGAL::to_double(e1 * e2)) > (1. - 1.e-5); +}; - auto G = G2.weld_vertices(); -#ifdef SVGFILL_DEBUG - obj << "o network_2\n"; - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - obj << "v " << CGAL::to_double(it->first.x()) << " " << CGAL::to_double(it->first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } - obj << std::flush; -#endif - - auto is_parallel_2degree_node = [](decltype(G)::vertex_const_iterator vit) { - auto it = vit->second.begin(); - auto& P = *it++; - auto& Q = *it++; - auto e1 = P - vit->first; - auto e2 = vit->first - Q; - if (e1.squared_length() == 0 || e2.squared_length() == 0) { - // @todo why does this happen? - return false; - } - e1 /= std::sqrt(CGAL::to_double(e1.squared_length())); - e2 /= std::sqrt(CGAL::to_double(e2.squared_length())); - return std::abs(CGAL::to_double(e1 * e2)) > (1. - 1.e-5); - }; - - { - // Remove colinear vertices - size_t n_vertices_removed = 0; - for (auto vit = G.vertices_begin(); vit != G.vertices_end();) { - if (vit->second.size() == 2) { - if (is_parallel_2degree_node(vit)) { - vit = G.eliminate_vertex(vit); - ++n_vertices_removed; - } else { - ++vit; - } +void eliminate_colinear_vertices(Graph2D& G) { + size_t n_vertices_removed = 0; + for (auto vit = G.vertices_begin(); vit != G.vertices_end();) { + if (vit->second.size() == 2) { + if (is_parallel_2degree_node(vit)) { + vit = G.eliminate_vertex(vit); + ++n_vertices_removed; } else { ++vit; } + } else { + ++vit; } - // std::cout << "Eliminated " << n_vertices_removed << " vertices" << std::endl; } +} - // Ortho edge slide - { - std::list> edges_to_remove, edges_to_insert; +void edge_slide(Graph2D& G) { + std::list> edges_to_remove, edges_to_insert; - for (auto vit = G.vertices_begin(); vit != G.vertices_end(); ++vit) { - auto& selected = vit->first; + for (auto vit = G.vertices_begin(); vit != G.vertices_end(); ++vit) { + auto& selected = vit->first; - if (vit->second.size() >= 3) { - for (auto vjt = vit->second.begin(); vjt != vit->second.end(); ++vjt) { - auto& neighbour = *vjt; - bool processed_neighbour = false; + if (vit->second.size() >= 3) { + for (auto vjt = vit->second.begin(); vjt != vit->second.end(); ++vjt) { + auto& neighbour = *vjt; + bool processed_neighbour = false; - if (G.find(neighbour)->second.size() == 2 && !is_parallel_2degree_node(G.find(neighbour))) { - auto vkt = G.find(neighbour)->second.begin(); - if (selected == *vkt) { - vkt++; - } - auto& other = *vkt; + if (G.find(neighbour)->second.size() == 2 && !is_parallel_2degree_node(G.find(neighbour))) { + auto vkt = G.find(neighbour)->second.begin(); + if (selected == *vkt) { + vkt++; + } + auto& other = *vkt; - if ((other - neighbour).squared_length() < (neighbour - selected).squared_length()) { - continue; - } + if ((other - neighbour).squared_length() < (neighbour - selected).squared_length()) { + continue; + } - auto incoming = CGAL::Ray_2(other, neighbour - other); - boost::optional> closest_neighbouring_segment; - boost::optional> closest_intersection_point; - K::FT sq_distance_along_ray = std::numeric_limits::infinity(); + auto incoming = CGAL::Ray_2(other, neighbour - other); + boost::optional> closest_neighbouring_segment; + boost::optional> closest_intersection_point; + K::FT sq_distance_along_ray = std::numeric_limits::infinity(); - for (auto vlt = vit->second.begin(); vlt != vit->second.end(); ++vlt) { - auto& other_neighbour = *vlt; - if (vlt != vjt) { - CGAL::Segment_2 neighbouring_segment(selected, other_neighbour); - auto x = CGAL::intersection(incoming, neighbouring_segment); - if (x) { - if (auto* xp = variant_get>(&*x)) { - auto dist = ((*xp) - other).squared_length(); - if (dist < sq_distance_along_ray) { - closest_neighbouring_segment = neighbouring_segment; - closest_intersection_point = *xp; - sq_distance_along_ray = dist; - } + for (auto vlt = vit->second.begin(); vlt != vit->second.end(); ++vlt) { + auto& other_neighbour = *vlt; + if (vlt != vjt) { + CGAL::Segment_2 neighbouring_segment(selected, other_neighbour); + auto x = CGAL::intersection(incoming, neighbouring_segment); + if (x) { + if (auto* xp = variant_get>(&*x)) { + auto dist = ((*xp) - other).squared_length(); + if (dist < sq_distance_along_ray) { + closest_neighbouring_segment = neighbouring_segment; + closest_intersection_point = *xp; + sq_distance_along_ray = dist; } } } } - - if (closest_intersection_point && closest_neighbouring_segment) { - edges_to_remove.push_back(*closest_neighbouring_segment); - edges_to_remove.push_back({ neighbour, selected }); - edges_to_insert.push_back({ closest_neighbouring_segment->source(), *closest_intersection_point }); - edges_to_insert.push_back({ closest_neighbouring_segment->target(), *closest_intersection_point }); - edges_to_insert.push_back({ neighbour, *closest_intersection_point }); - - processed_neighbour = true; - } } - if (processed_neighbour) { - // Only one neigbour is processed because otherwise we obtain intersections - break; + + if (closest_intersection_point && closest_neighbouring_segment) { + edges_to_remove.push_back(*closest_neighbouring_segment); + edges_to_remove.push_back({neighbour, selected}); + edges_to_insert.push_back({closest_neighbouring_segment->source(), *closest_intersection_point}); + edges_to_insert.push_back({closest_neighbouring_segment->target(), *closest_intersection_point}); + edges_to_insert.push_back({neighbour, *closest_intersection_point}); + + processed_neighbour = true; } } + if (processed_neighbour) { + // Only one neigbour is processed because otherwise we obtain intersections + break; + } } } - - for (auto& s : edges_to_remove) { - G.remove_edge(s.source(), s.target()); - } - - - for (auto& s : edges_to_insert) { - G.insert(s.source(), s.target()); - } - -#ifdef SVGFILL_DEBUG - obj << "o network_3\n"; - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - obj << "v " << CGAL::to_double(it->first.x()) << " " << CGAL::to_double(it->first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } -#endif } - // Now plot the edges on an arrangement in order to find planar cycles - // and merge the corridor-halves with their neighbouring input polygon - - Arrangement_2 arr; - - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - if (it->first == it->second) { - continue; - } - CGAL::insert(arr, Segment_2(it->first, it->second)); + for (auto& s : edges_to_remove) { + G.remove_edge(s.source(), s.target()); } - std::list> move_ops; - std::list> edge_ops; + for (auto& s : edges_to_insert) { + G.insert(s.source(), s.target()); + } +} + +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 +){ + std::list> constructed_segments; for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) { if (it->second.size() == 1) { auto& M = it->first; - decltype(midpoint_to_segment)::mapped_type* q = nullptr; + const std::pair* q = nullptr; if (midpoint_to_segment.find(M) == midpoint_to_segment.end()) { typename K::FT min_sq_distance = std::numeric_limits::infinity(); @@ -1299,7 +1030,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } else { - q = &midpoint_to_segment[M]; + q = &midpoint_to_segment.find(M)->second; } if (q == nullptr) { @@ -1309,7 +1040,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v bool handled_as_graph_path = false; // distance from unioned - shoot ray? - if (segment_to_input_facet[*q].size() == 2) { + 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_bounded_side(M)) { @@ -1339,10 +1070,10 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v Graph2D GGG(bnd); GGG.refine(*GGG.query(*closest_intersection_point, 0.01), *closest_intersection_point); - std::array>, 2> input_points = { { {}, {} } }; + std::array>, 2> input_points = {{{}, {}}}; size_t i = 0; - for (auto& fac : segment_to_input_facet[*q]) { + for (auto& fac : segment_to_input_facet.find(*q)->second) { for (auto it = fac->vertices_begin(); it != fac->vertices_end(); ++it) { auto seg = GGG.query(*it, 0.01); if (seg) { @@ -1361,13 +1092,13 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (!a1.empty() && !a2.empty()) { if (M != *closest_intersection_point) { - edge_ops.push_front({ M, *closest_intersection_point }); + constructed_segments.push_front({M, *closest_intersection_point}); } for (auto it = a1.begin(); it != a1.end() && std::next(it) != a1.end(); ++it) { - edge_ops.push_front({ *it, *(std::next(it)) }); + constructed_segments.push_front({*it, *(std::next(it))}); } for (auto it = a2.begin(); it != a2.end() && std::next(it) != a2.end(); ++it) { - edge_ops.push_front({ *it, *(std::next(it)) }); + constructed_segments.push_front({*it, *(std::next(it))}); } handled_as_graph_path = true; @@ -1381,8 +1112,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (!handled_as_graph_path) { // else we choose to map point to the midpoint of the found two close points. - auto pq = close_input_point(q->first); - auto pr = close_input_point(q->second); + auto pq = segment_lookup.close_input_point(q->first); + auto pr = segment_lookup.close_input_point(q->second); auto Q = pq.second; auto R = pr.second; @@ -1392,16 +1123,16 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // where Q and R are co-located, because the point R' is further away // in that case M + M-Q should gives is x that we then project onto the // input boundary - // - // - // ┌───────┐ - // │ │ - // │ │ - // │ │ - // └───────o <--Q,R - // - // ────────o <--M - // + // + // + // ┌───────┐ + // │ │ + // │ │ + // │ │ + // └───────o <--Q,R + // + // ────────o <--M + // // ┌───────x───────────────o <---R' // │ │ // │ │ @@ -1410,93 +1141,22 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // └───────────────────────┘ // @todo is this projection actually necessary or is it already 'exact enough'? - R = project_input_point(M + (M - Q)).second; + R = segment_lookup.project_input_point(M + (M - Q)).second; } auto avg = CGAL::ORIGIN + ((Q - CGAL::ORIGIN) + (R - CGAL::ORIGIN)) / 2; - move_ops.push_front({ M, avg }); - edge_ops.push_front({ avg, Q }); - edge_ops.push_front({ avg, R }); - + constructed_segments.push_front({M, avg}); + constructed_segments.push_front({avg, Q}); + constructed_segments.push_front({avg, R}); } } } -#ifdef SVGFILL_DEBUG - obj << "o network_4\n"; -#endif + return constructed_segments; +} - // note that we actually don't move but draw an edge - for (auto& pq : move_ops) { - if (pq.first == pq.second) { - continue; - } - CGAL::insert(arr, Segment_2(pq.first, pq.second)); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(pq.first.x()) << " " << CGAL::to_double(pq.first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(pq.second.x()) << " " << CGAL::to_double(pq.second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; -#endif - } - - - for (auto& pq : edge_ops) { - if (pq.first == pq.second) { - continue; - } - CGAL::insert(arr, Segment_2(pq.first, pq.second)); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(pq.first.x()) << " " << CGAL::to_double(pq.first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(pq.second.x()) << " " << CGAL::to_double(pq.second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; -#endif - } - - // Plot input polygons - 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))); - } - } - - -#ifdef SVGFILL_DEBUG - { - obj << "o arrangement_1\n"; - for (auto it = arr.edges_begin(); it != arr.edges_end(); ++it) { - auto& p = it->source()->point(); - auto& q = it->target()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - obj << "v " << CGAL::to_double(q.x()) << " " << CGAL::to_double(q.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } - } -#endif - - /* { - // debug, add outer bounds so that we can plot the face for any remaining edges - auto poly = unioned_polygons.front().outer_boundary(); - for (size_t i = 0; i != poly.size(); ++i) { - auto j = (i + 1) % poly.size(); - CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); - } - } */ - - // Now loop over the arrangement faces, when a face coincides with a point on the - // 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. - - size_t face_id = 0; +void fuse_corridor_halves_with_input(Arrangement_2& arr, Graph2D& G, SegmentLookup& segment_lookup, const Polygon_list& input_polygons, DebugWriter& debug_output) { std::set edges_to_remove; for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { @@ -1535,7 +1195,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto& p = curr->source()->point(); auto& q = curr->target()->point(); auto center = CGAL::ORIGIN + (((p - CGAL::ORIGIN) + (q - CGAL::ORIGIN)) / 2); - auto p1index = input_polygon_boundary(center); + auto p1index = segment_lookup.input_polygon_boundary(center); const bool on_orig_bound = p1index != input_polygons.end(); if (on_orig_bound) { if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { @@ -1553,7 +1213,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto& p = curr->source()->point(); auto& q = curr->target()->point(); auto center = CGAL::ORIGIN + (((p - CGAL::ORIGIN) + (q - CGAL::ORIGIN)) / 2); - auto p1index = input_polygon_boundary(center); + auto p1index = segment_lookup.input_polygon_boundary(center); const bool on_orig_bound = p1index != input_polygons.end(); if (on_orig_bound) { if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { @@ -1569,126 +1229,314 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } - -#ifdef SVGFILL_DEBUG - write_polygon_to_svg(svg, circ_to_poly(it->outer_ccb())); - - obj << "o " << "face_"; - if (is_corridor) { - obj << "corri_"; - } - obj << face_id++ << "\n"; - - std::ostringstream oss; - - { - auto vv = vi; - auto curr = it->outer_ccb(); - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == it->outer_ccb()) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != it->outer_ccb()); - } - - for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { - auto vv = vi; - auto curr = *jt; - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == *jt) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != *jt); - } - - obj << oss.str(); -#endif } size_t remove_id = 0; for (auto& e : edges_to_remove) { -#ifdef SVGFILL_DEBUG - obj << "o " << "remove_" << remove_id++ << "\n"; - { - auto& p = e->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - } - { - auto& p = e->target()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - } - obj << "l " << vi++; - obj << " " << vi++ << std::endl; -#endif + debug_output.write_segment(e->source()->point(), e->target()->point(), "arr_remove_edge_" + std::to_string(remove_id++)); CGAL::remove_edge(arr, e); } +} + +class timer { + class entry { + public: + 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; + } + + private: + std::map::const_iterator start_it; + }; + + public: + entry start(const std::string& name) { + return timings_.insert({name, std::chrono::high_resolution_clock::now()}).first; + } + + private: + std::map< + std::string, + std::chrono::high_resolution_clock::time_point> + timings_; +}; + +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; + // 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); + + 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 + + timer timer; + + auto t0 = timer.start("input"); + + debug_output.write_polygons(input_polygons_, "input"); + + if (polygon_offset_distance < 0.) { + polygon_offset_distance = estimate_polygon_offset_distance(input_polygons_); + } + + // Create copy to make mutable for cleaning + auto input_polygons = input_polygons_; + + for (auto& polygon : input_polygons) { + clean_polygon(polygon); + } + + { + decltype(input_polygons) split_polygons; + for (auto& poly : input_polygons) { + split_self_intersecting_polygon(poly, std::back_inserter(split_polygons)); + } + std::swap(input_polygons, split_polygons); + } + + t0.stop(); + t0 = timer.start("overlap elimination"); + + eliminate_overlaps(OVERLAP_RESOLUTION_DISTANCE, input_polygons); + + t0.stop(); + + // [NB Nov 6] we cannot do this anymore because it could revert the spacing between input polygons + // that touch in the corner. + // Now that overlaps/touches at corners are handled more locally only a small indent is produced + // which would be undone by means of an inset+offset. + // + // [NB Nov 10] this is actually still necessary though, but we apply a much smaller distance now + // to keep the overlap eliminations in tact + // + // 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); + } + + debug_output.write_polygons(input_polygons, "processed_input"); + + SegmentLookup segment_lookup(input_polygons); + + 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 inner_offset = 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(inner_offset, "outer_perimiter"); + + 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) { + std::vector working_copy; + working_copy.emplace_back(i); + + for (auto& r : input_polygons) { + std::vector temp_working_copy; + for (auto& wc : working_copy) { + CGAL::difference(wc, r, std::back_inserter(temp_working_copy)); + } + working_copy = temp_working_copy; + } + difference_result.insert(difference_result.end(), working_copy.begin(), working_copy.end()); + } + + t0.stop(); + t0 = timer.start("corridor triangulation"); + + // subdivide difference_result to have better more detailed triangulation and therefore less-pronounced artefacts in midpoint network + + for (auto& pwh : difference_result) { + difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 8., pwh)); + // difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 64., pwh)); + } + + debug_output.write_polygons(difference_result_subdivided, "corridor_subdivided"); + + std::vector> triangular_polygons; + for (auto& pwh : difference_result_subdivided) { + CGAL::Polygon_triangulation_decomposition_2 decompositor; + decompositor(pwh, std::back_inserter(triangular_polygons)); + } + + t0.stop(); + + /* + * // @todo decide whether this is smart or not + * // Would this not hurt topology too much? + triangular_polygons.erase(std::remove_if(triangular_polygons.begin(), triangular_polygons.end(), [](const CGAL::Polygon_2& p) { + return CGAL::to_double(p.area()) < 1.e-8; + }), triangular_polygons.end()); + */ + + t0 = timer.start("center line"); + + debug_output.write_polygons(triangular_polygons, "triangulated_corridor"); + + 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) { + debug_output.write_segment(p.first, q, "network_1"); + } + } + + 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 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(); + + 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); + + edge_slide(G); + + 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, inner_offset, 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 + Arrangement_2 arr; + G.to_arrangement(arr); + + for (auto& pq : segments) { + if (pq.first == pq.second) { + continue; + } + CGAL::insert(arr, Segment_2(pq.first, pq.second)); + + debug_output.write_segment(pq.first, pq.second, "extended_segments"); + } + + // Write input polygons to arrangement_2 + 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))); + } + } + + // Just for the automatic numbering, create a full vector + std::vector temp; + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (it->is_unbounded()) { + continue; + } + temp.push_back(circ_to_poly(it->outer_ccb())); + } + debug_output.write_polygons(temp, "arr_faces"); + + + /* { + // debug, add outer bounds so that we can plot the face for any remaining edges + auto poly = unioned_polygons.front().outer_boundary(); + for (size_t i = 0; i != poly.size(); ++i) { + auto j = (i + 1) % poly.size(); + CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); + } + } */ + + // Now loop over the arrangement faces, when a face coincides with a point on the + // 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. + + fuse_corridor_halves_with_input(arr, G, segment_lookup, input_polygons, debug_output); + + t0.stop(); for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { if (it->is_unbounded()) { continue; } - output_polygons.push_back(circ_to_poly(it->outer_ccb())); - -#ifdef SVGFILL_DEBUG - write_polygon_to_svg(svg, circ_to_poly(it->outer_ccb())); - - obj << "o " << "merged_face_"; - obj << face_id++ << "\n"; - - std::ostringstream oss; - - { - auto vv = vi; - auto curr = it->outer_ccb(); - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == it->outer_ccb()) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != it->outer_ccb()); - } - - for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { - auto vv = vi; - auto curr = *jt; - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == *jt) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != *jt); - } - - obj << oss.str(); -#endif } -#ifdef SVGFILL_DEBUG - svg << "\n"; -#endif + debug_output.write_polygons(output_polygons, "arr_faces_merged"); } #ifndef SVGFILL_MAIN diff --git a/src/svgfill/src/graph_2d.h b/src/svgfill/src/graph_2d.h index aaaa059b4e..da2b5014ec 100644 --- a/src/svgfill/src/graph_2d.h +++ b/src/svgfill/src/graph_2d.h @@ -2,8 +2,10 @@ #define GRAPH_2D_H #ifdef SVGFILL_DEBUG +#if 0 #include #endif +#endif template class Graph2D { @@ -334,6 +336,16 @@ public: return Graph2D(input_adjacency_list); } + template + void to_arrangement(T& arr) { + for (auto it = edges_begin(); it != edges_end(); ++it) { + if (it->first == it->second) { + continue; + } + CGAL::insert(arr, CGAL::Segment_2(it->first, it->second)); + } + } + void assert_symmetric() { #ifdef SVGFILL_DEBUG #if 0