From 246ff2653dcc26d042c8205a26c6c22b36b768ab Mon Sep 17 00:00:00 2001 From: Arthur Date: Fri, 10 Mar 2023 20:44:07 +0800 Subject: [PATCH] ENH: improve bridge direction detection Thanks to Prusa. Jira: STUDIO-2453 Change-Id: Iadcc59df44d5abc552f5d558a500fd9bcd66d43f (cherry picked from commit c19156fd037df4231f3e0cb1e9a899c9b7525372) --- src/libslic3r/BridgeDetector.cpp | 3 +- src/libslic3r/BridgeDetector.hpp | 61 ++++++++++- src/libslic3r/CMakeLists.txt | 2 + src/libslic3r/LayerRegion.cpp | 10 ++ src/libslic3r/PrincipalComponents2D.cpp | 140 ++++++++++++++++++++++++ src/libslic3r/PrincipalComponents2D.hpp | 24 ++++ 6 files changed, 238 insertions(+), 2 deletions(-) create mode 100644 src/libslic3r/PrincipalComponents2D.cpp create mode 100644 src/libslic3r/PrincipalComponents2D.hpp diff --git a/src/libslic3r/BridgeDetector.cpp b/src/libslic3r/BridgeDetector.cpp index a60f45f11..09a2c2dc3 100644 --- a/src/libslic3r/BridgeDetector.cpp +++ b/src/libslic3r/BridgeDetector.cpp @@ -160,7 +160,8 @@ bool BridgeDetector::detect_angle(double bridge_direction_override) // if any other direction is within extrusion width of coverage, prefer it if shorter // TODO: There are two options here - within width of the angle with most coverage, or within width of the currently perferred? size_t i_best = 0; - for (size_t i = 1; i < candidates.size() && abs(candidates[i_best].archored_percent - candidates[i].archored_percent) < EPSILON; ++ i) +// for (size_t i = 1; i < candidates.size() && abs(candidates[i_best].archored_percent - candidates[i].archored_percent) < EPSILON; ++ i) + for (size_t i = 1; i < candidates.size() && candidates[i_best].coverage - candidates[i].coverage < this->spacing; ++ i) if (candidates[i].max_length < candidates[i_best].max_length) i_best = i; diff --git a/src/libslic3r/BridgeDetector.hpp b/src/libslic3r/BridgeDetector.hpp index d422949cd..064e40dd5 100644 --- a/src/libslic3r/BridgeDetector.hpp +++ b/src/libslic3r/BridgeDetector.hpp @@ -1,6 +1,12 @@ #ifndef slic3r_BridgeDetector_hpp_ #define slic3r_BridgeDetector_hpp_ +#include "ClipperUtils.hpp" +#include "Line.hpp" +#include "Point.hpp" +#include "Polygon.hpp" +#include "Polyline.hpp" +#include "PrincipalComponents2D.hpp" #include "libslic3r.h" #include "ExPolygon.hpp" #include @@ -48,7 +54,7 @@ private: // the best direction is the one causing most lines to be bridged (thus most coverage) bool operator<(const BridgeDirection &other) const { // Initial sort by coverage only - comparator must obey strict weak ordering - return this->archored_percent > other.archored_percent; + return this->coverage > other.coverage;//this->archored_percent > other.archored_percent; }; double angle; double coverage; @@ -65,6 +71,59 @@ private: ExPolygons _anchor_regions; }; + +//return ideal bridge direction and unsupported bridge endpoints distance. +inline std::tuple detect_bridging_direction(const Polygons &to_cover, const Polygons &anchors_area) +{ + Polygons overhang_area = diff(to_cover, anchors_area); + Polylines floating_polylines = diff_pl(to_polylines(overhang_area), expand(anchors_area, float(SCALED_EPSILON))); + + if (floating_polylines.empty()) { + // consider this area anchored from all sides, pick bridging direction that will likely yield shortest bridges + auto [pc1, pc2] = compute_principal_components(overhang_area); + if (pc2 == Vec2f::Zero()) { // overhang may be smaller than resolution. In this case, any direction is ok + return {Vec2d{1.0,0.0}, 0.0}; + } else { + return {pc2.normalized().cast(), 0.0}; + } + } + + // Overhang is not fully surrounded by anchors, in that case, find such direction that will minimize the number of bridge ends/180turns in the air + Lines floating_edges = to_lines(floating_polylines); + std::unordered_map directions{}; + for (const Line &l : floating_edges) { + Vec2d normal = l.normal().cast().normalized(); + double quantized_angle = std::ceil(std::atan2(normal.y(),normal.x()) * 1000.0); + directions.emplace(quantized_angle, normal); + } + std::vector> direction_costs{}; + // it is acutally cost of a perpendicular bridge direction - we find the minimal cost and then return the perpendicular dir + for (const auto& d : directions) { + direction_costs.emplace_back(d.second, 0.0); + } + + for (const Line &l : floating_edges) { + Vec2d line = (l.b - l.a).cast(); + for (auto &dir_cost : direction_costs) { + // the dot product already contains the length of the line. dir_cost.first is normalized. + dir_cost.second += std::abs(line.dot(dir_cost.first)); + } + } + + Vec2d result_dir = Vec2d::Ones(); + double min_cost = std::numeric_limits::max(); + for (const auto &cost : direction_costs) { + if (cost.second < min_cost) { + // now flip the orientation back and return the direction of the bridge extrusions + result_dir = Vec2d{cost.first.y(), -cost.first.x()}; + min_cost = cost.second; + } + } + + return {result_dir, min_cost}; +}; + + } #endif diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 43d10b5a2..419f1c245 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -220,6 +220,8 @@ set(lisbslic3r_sources PresetBundle.hpp ProjectTask.cpp ProjectTask.hpp + PrincipalComponents2D.hpp + PrincipalComponents2D.cpp AppConfig.cpp AppConfig.hpp Print.cpp diff --git a/src/libslic3r/LayerRegion.cpp b/src/libslic3r/LayerRegion.cpp index f66ba002c..63e32a8a9 100644 --- a/src/libslic3r/LayerRegion.cpp +++ b/src/libslic3r/LayerRegion.cpp @@ -305,6 +305,15 @@ void LayerRegion::process_external_surfaces(const Layer *lower_layer, const Poly // would get merged into a single one while they need different directions // also, supply the original expolygon instead of the grown one, because in case // of very thin (but still working) anchors, the grown expolygon would go beyond them + double custom_angle = Geometry::deg2rad(this->region().config().bridge_angle.value); + if (custom_angle > 0.0) { + bridges[idx_last].bridge_angle = custom_angle; + } else { + auto [bridging_dir, unsupported_dist] = detect_bridging_direction(to_polygons(initial), to_polygons(lower_layer->lslices)); + bridges[idx_last].bridge_angle = PI + std::atan2(bridging_dir.y(), bridging_dir.x()); + } + + /* BridgeDetector bd(initial, lower_layer->lslices, this->bridging_flow(frInfill, object_config.thick_bridges).scaled_width()); #ifdef SLIC3R_DEBUG printf("Processing bridge at layer %zu:\n", this->layer()->id()); @@ -321,6 +330,7 @@ void LayerRegion::process_external_surfaces(const Layer *lower_layer, const Poly // using a bridging flow, therefore it makes sense to respect the custom bridging direction. bridges[idx_last].bridge_angle = custom_angle; } + */ // without safety offset, artifacts are generated (GH #2494) surfaces_append(bottom, union_safety_offset_ex(grown), bridges[idx_last]); } diff --git a/src/libslic3r/PrincipalComponents2D.cpp b/src/libslic3r/PrincipalComponents2D.cpp new file mode 100644 index 000000000..7bdf79315 --- /dev/null +++ b/src/libslic3r/PrincipalComponents2D.cpp @@ -0,0 +1,140 @@ +#include "PrincipalComponents2D.hpp" +#include "Point.hpp" + +namespace Slic3r { + + + +// returns triangle area, first_moment_of_area_xy, second_moment_of_area_xy, second_moment_of_area_covariance +// none of the values is divided/normalized by area. +// The function computes intgeral over the area of the triangle, with function f(x,y) = x for first moments of area (y is analogous) +// f(x,y) = x^2 for second moment of area +// and f(x,y) = x*y for second moment of area covariance +std::tuple compute_moments_of_area_of_triangle(const Vec2f &a, const Vec2f &b, const Vec2f &c) +{ + // based on the following guide: + // Denote the vertices of S by a, b, c. Then the map + // g:(u,v)↦a+u(b−a)+v(c−a) , + // which in coordinates appears as + // g:(u,v)↦{x(u,v)y(u,v)=a1+u(b1−a1)+v(c1−a1)=a2+u(b2−a2)+v(c2−a2) ,(1) + // obviously maps S′ bijectively onto S. Therefore the transformation formula for multiple integrals steps into action, and we obtain + // ∫Sf(x,y)d(x,y)=∫S′f(x(u,v),y(u,v))∣∣Jg(u,v)∣∣ d(u,v) . + // In the case at hand the Jacobian determinant is a constant: From (1) we obtain + // Jg(u,v)=det[xuyuxvyv]=(b1−a1)(c2−a2)−(c1−a1)(b2−a2) . + // Therefore we can write + // ∫Sf(x,y)d(x,y)=∣∣Jg∣∣∫10∫1−u0f~(u,v) dv du , + // where f~ denotes the pullback of f to S′: + // f~(u,v):=f(x(u,v),y(u,v)) . + // Don't forget taking the absolute value of Jg! + + float jacobian_determinant_abs = std::abs((b.x() - a.x()) * (c.y() - a.y()) - (c.x() - a.x()) * (b.y() - a.y())); + + // coordinate transform: gx(u,v) = a.x + u * (b.x - a.x) + v * (c.x - a.x) + // coordinate transform: gy(u,v) = a.y + u * (b.y - a.y) + v * (c.y - a.y) + // second moment of area for x: f(x, y) = x^2; + // f(gx(u,v), gy(u,v)) = gx(u,v)^2 = ... (long expanded form) + + // result is Int_T func = jacobian_determinant_abs * Int_0^1 Int_0^1-u func(gx(u,v), gy(u,v)) dv du + // integral_0^1 integral_0^(1 - u) (a + u (b - a) + v (c - a))^2 dv du = 1/12 (a^2 + a (b + c) + b^2 + b c + c^2) + + Vec2f second_moment_of_area_xy = jacobian_determinant_abs * + (a.cwiseProduct(a) + b.cwiseProduct(b) + b.cwiseProduct(c) + c.cwiseProduct(c) + + a.cwiseProduct(b + c)) / + 12.0f; + // second moment of area covariance : f(x, y) = x*y; + // f(gx(u,v), gy(u,v)) = gx(u,v)*gy(u,v) = ... (long expanded form) + //(a_1 + u * (b_1 - a_1) + v * (c_1 - a_1)) * (a_2 + u * (b_2 - a_2) + v * (c_2 - a_2)) + // == (a_1 + u (b_1 - a_1) + v (c_1 - a_1)) (a_2 + u (b_2 - a_2) + v (c_2 - a_2)) + + // intermediate result: integral_0^(1 - u) (a_1 + u (b_1 - a_1) + v (c_1 - a_1)) (a_2 + u (b_2 - a_2) + v (c_2 - a_2)) dv = + // 1/6 (u - 1) (-c_1 (u - 1) (a_2 (u - 1) - 3 b_2 u) - c_2 (u - 1) (a_1 (u - 1) - 3 b_1 u + 2 c_1 (u - 1)) + 3 b_1 u (a_2 (u - 1) - 2 + // b_2 u) + a_1 (u - 1) (3 b_2 u - 2 a_2 (u - 1))) result = integral_0^1 1/6 (u - 1) (-c_1 (u - 1) (a_2 (u - 1) - 3 b_2 u) - c_2 (u - + // 1) (a_1 (u - 1) - 3 b_1 u + 2 c_1 (u - 1)) + 3 b_1 u (a_2 (u - 1) - 2 b_2 u) + a_1 (u - 1) (3 b_2 u - 2 a_2 (u - 1))) du = + // 1/24 (a_2 (b_1 + c_1) + a_1 (2 a_2 + b_2 + c_2) + b_2 c_1 + b_1 c_2 + 2 b_1 b_2 + 2 c_1 c_2) + // result is Int_T func = jacobian_determinant_abs * Int_0^1 Int_0^1-u func(gx(u,v), gy(u,v)) dv du + float second_moment_of_area_covariance = jacobian_determinant_abs * (1.0f / 24.0f) * + (a.y() * (b.x() + c.x()) + a.x() * (2.0f * a.y() + b.y() + c.y()) + b.y() * c.x() + + b.x() * c.y() + 2.0f * b.x() * b.y() + 2.0f * c.x() * c.y()); + + float area = jacobian_determinant_abs * 0.5f; + + Vec2f first_moment_of_area_xy = jacobian_determinant_abs * (a + b + c) / 6.0f; + + return {area, first_moment_of_area_xy, second_moment_of_area_xy, second_moment_of_area_covariance}; +}; + +// returns two eigenvectors of the area covered by given polygons. The vectors are sorted by their corresponding eigenvalue, largest first +std::tuple compute_principal_components(const Polygons &polys) +{ + Vec2f centroid_accumulator = Vec2f::Zero(); + Vec2f second_moment_of_area_accumulator = Vec2f::Zero(); + float second_moment_of_area_covariance_accumulator = 0.0f; + float area = 0.0f; + + for (const Polygon &poly : polys) { + Vec2f p0 = unscaled(poly.first_point()).cast(); + for (size_t i = 2; i < poly.points.size(); i++) { + Vec2f p1 = unscaled(poly.points[i - 1]).cast(); + Vec2f p2 = unscaled(poly.points[i]).cast(); + + float sign = cross2(p1 - p0, p2 - p1) > 0 ? 1.0f : -1.0f; + + auto [triangle_area, first_moment_of_area, second_moment_area, + second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2); + area += sign * triangle_area; + centroid_accumulator += sign * first_moment_of_area; + second_moment_of_area_accumulator += sign * second_moment_area; + second_moment_of_area_covariance_accumulator += sign * second_moment_of_area_covariance; + } + } + + if (area <= 0.0) { + return {Vec2f::Zero(), Vec2f::Zero()}; + } + + Vec2f centroid = centroid_accumulator / area; + Vec2f variance = second_moment_of_area_accumulator / area - centroid.cwiseProduct(centroid); + double covariance = second_moment_of_area_covariance_accumulator / area - centroid.x() * centroid.y(); +#if 0 + std::cout << "area : " << area << std::endl; + std::cout << "variancex : " << variance.x() << std::endl; + std::cout << "variancey : " << variance.y() << std::endl; + std::cout << "covariance : " << covariance << std::endl; +#endif + if (abs(covariance) < EPSILON) { + std::tuple result{Vec2f{variance.x(), 0.0}, Vec2f{0.0, variance.y()}}; + if (variance.y() > variance.x()) { + return {std::get<1>(result), std::get<0>(result)}; + } else + return result; + } + + // now we find the first principal component of the covered area by computing max eigenvalue and the correspoding eigenvector of + // covariance matrix + // covaraince matrix C is : | VarX Cov | + // | Cov VarY | + // Eigenvalues are solutions to det(C - lI) = 0, where l is the eigenvalue and I unit matrix + // Eigenvector for eigenvalue l is any vector v such that Cv = lv + + float eigenvalue_a = 0.5f * (variance.x() + variance.y() + + sqrt((variance.x() - variance.y()) * (variance.x() - variance.y()) + 4.0f * covariance * covariance)); + float eigenvalue_b = 0.5f * (variance.x() + variance.y() - + sqrt((variance.x() - variance.y()) * (variance.x() - variance.y()) + 4.0f * covariance * covariance)); + Vec2f eigenvector_a{(eigenvalue_a - variance.y()) / covariance, 1.0f}; + Vec2f eigenvector_b{(eigenvalue_b - variance.y()) / covariance, 1.0f}; + +#if 0 + std::cout << "eigenvalue_a: " << eigenvalue_a << std::endl; + std::cout << "eigenvalue_b: " << eigenvalue_b << std::endl; + std::cout << "eigenvectorA: " << eigenvector_a.x() << " " << eigenvector_a.y() << std::endl; + std::cout << "eigenvectorB: " << eigenvector_b.x() << " " << eigenvector_b.y() << std::endl; +#endif + + if (eigenvalue_a > eigenvalue_b) { + return {eigenvector_a, eigenvector_b}; + } else { + return {eigenvector_b, eigenvector_a}; + } +} + +} \ No newline at end of file diff --git a/src/libslic3r/PrincipalComponents2D.hpp b/src/libslic3r/PrincipalComponents2D.hpp new file mode 100644 index 000000000..dc8897a7a --- /dev/null +++ b/src/libslic3r/PrincipalComponents2D.hpp @@ -0,0 +1,24 @@ +#ifndef slic3r_PrincipalComponents2D_hpp_ +#define slic3r_PrincipalComponents2D_hpp_ + +#include "AABBTreeLines.hpp" +#include "BoundingBox.hpp" +#include "libslic3r.h" +#include +#include "Polygon.hpp" + +namespace Slic3r { + +// returns triangle area, first_moment_of_area_xy, second_moment_of_area_xy, second_moment_of_area_covariance +// none of the values is divided/normalized by area. +// The function computes intgeral over the area of the triangle, with function f(x,y) = x for first moments of area (y is analogous) +// f(x,y) = x^2 for second moment of area +// and f(x,y) = x*y for second moment of area covariance +std::tuple compute_moments_of_area_of_triangle(const Vec2f &a, const Vec2f &b, const Vec2f &c); + +// returns two eigenvectors of the area covered by given polygons. The vectors are sorted by their corresponding eigenvalue, largest first +std::tuple compute_principal_components(const Polygons &polys); + +} + +#endif \ No newline at end of file