From 782fa47791080a53336c3af0fddc6a51eb282c8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20=C5=A0ach?= Date: Wed, 4 Oct 2023 13:46:12 +0200 Subject: [PATCH 1/2] Remove usage of polygons_covered_by_width. In SupportSpotsGenerator there is a need to compute integral over extrusions. The mentioned function is slow for this use case. As it solves a harder problem. It is better to iterate the extrusion lines directly. --- src/libslic3r/SupportSpotsGenerator.cpp | 71 ++++++++++++++++--- src/libslic3r/SupportSpotsGenerator.hpp | 4 ++ .../test_support_spots_generator.cpp | 15 +++- 3 files changed, 79 insertions(+), 11 deletions(-) diff --git a/src/libslic3r/SupportSpotsGenerator.cpp b/src/libslic3r/SupportSpotsGenerator.cpp index 679ef0215f..7ca54793ec 100644 --- a/src/libslic3r/SupportSpotsGenerator.cpp +++ b/src/libslic3r/SupportSpotsGenerator.cpp @@ -174,6 +174,32 @@ Integrals::Integrals (const Polygons& polygons) { } } +Integrals::Integrals(const Polylines& polylines, const std::vector& widths) { + assert(extrusion_lines.size() == widths.size()); + for (size_t i = 0; i < polylines.size(); ++i) { + Lines polyline{polylines[i].lines()}; + float width{widths[i]}; + for (const Line& line : polyline) { + Vec2f line_direction = unscaled(line.vector()).cast(); + Vec2f normal{line_direction.y(), -line_direction.x()}; + normal.normalize(); + + Vec2f line_a = unscaled(line.a).cast(); + Vec2f line_b = unscaled(line.b).cast(); + Vec2crd a = scaled(Vec2f{line_a + normal * width/2}); + Vec2crd b = scaled(Vec2f{line_b + normal * width/2}); + Vec2crd c = scaled(Vec2f{line_b - normal * width/2}); + Vec2crd d = scaled(Vec2f{line_a - normal * width/2}); + + const Polygon ractangle({a, b, c, d}); + Integrals integrals{{ractangle}}; + this->area += integrals.area; + this->x_i += integrals.x_i; + this->x_i_squared += integrals.x_i_squared; + this->xy += integrals.xy; + } + } +} SliceConnection estimate_slice_connection(size_t slice_idx, const Layer *layer) { @@ -473,18 +499,43 @@ ObjectPart::ObjectPart( continue; } - const Polygons polygons{collection->polygons_covered_by_width()}; + for (const ExtrusionEntity* entity: collection->flatten()) { + Polylines polylines; + std::vector widths; - const Integrals integrals{polygons}; - const float volume = integrals.area * layer_height; - this->volume += volume; - this->volume_centroid_accumulator += to_3d(integrals.x_i, center_z * integrals.area) / integrals.area * volume; + const auto* path = dynamic_cast(entity); + if (path !=nullptr) { + polylines.push_back(path->as_polyline()); + widths.push_back(path->width()); + } - if (this->connected_to_bed) { - this->sticking_area += integrals.area; - this->sticking_centroid_accumulator += to_3d(integrals.x_i, bottom_z * integrals.area); - this->sticking_second_moment_of_area_accumulator += integrals.x_i_squared; - this->sticking_second_moment_of_area_covariance_accumulator += integrals.xy; + const auto* loop = dynamic_cast(entity); + if (loop !=nullptr) { + for (const ExtrusionPath& path : loop->paths) { + polylines.push_back(path.as_polyline()); + widths.push_back(path.width()); + } + } + + const auto* multi_path = dynamic_cast(entity); + if (multi_path !=nullptr) { + for (const ExtrusionPath& path : multi_path->paths) { + polylines.push_back(path.as_polyline()); + widths.push_back(path.width()); + } + } + + const Integrals integrals{polylines, widths}; + const float volume = integrals.area * layer_height; + this->volume += volume; + this->volume_centroid_accumulator += to_3d(integrals.x_i, center_z * integrals.area) / integrals.area * volume; + + if (this->connected_to_bed) { + this->sticking_area += integrals.area; + this->sticking_centroid_accumulator += to_3d(integrals.x_i, bottom_z * integrals.area); + this->sticking_second_moment_of_area_accumulator += integrals.x_i_squared; + this->sticking_second_moment_of_area_covariance_accumulator += integrals.xy; + } } } diff --git a/src/libslic3r/SupportSpotsGenerator.hpp b/src/libslic3r/SupportSpotsGenerator.hpp index f13d0c9bb4..0f08550b7a 100644 --- a/src/libslic3r/SupportSpotsGenerator.hpp +++ b/src/libslic3r/SupportSpotsGenerator.hpp @@ -151,6 +151,7 @@ class Integrals{ * @param polygons List of polygons specifing the domain. */ explicit Integrals(const Polygons& polygons); + explicit Integrals(const Polylines& polylines, const std::vector& widths); // TODO refactor and delete the default constructor Integrals() = default; @@ -159,6 +160,9 @@ class Integrals{ Vec2f x_i{Vec2f::Zero()}; Vec2f x_i_squared{Vec2f::Zero()}; float xy{}; + +private: + void add(const Integrals& other); }; float compute_second_moment( diff --git a/tests/libslic3r/test_support_spots_generator.cpp b/tests/libslic3r/test_support_spots_generator.cpp index 0cd8c51f8e..1b51ccb0ac 100644 --- a/tests/libslic3r/test_support_spots_generator.cpp +++ b/tests/libslic3r/test_support_spots_generator.cpp @@ -6,7 +6,7 @@ using namespace Slic3r; using namespace SupportSpotsGenerator; -TEST_CASE("Numerical integral calculation compared with exact solution.", "[SupportSpotsGenerator]") { +TEST_CASE("Numerical integral over polygon calculation compared with exact solution.", "[SupportSpotsGenerator]") { const float width = 10; const float height = 20; const Polygon polygon = { @@ -24,6 +24,19 @@ TEST_CASE("Numerical integral calculation compared with exact solution.", "[Supp CHECK(integrals.x_i_squared.y() == Approx(width * std::pow(height, 3) / 12)); } +TEST_CASE("Numerical integral over line calculation compared with exact solution.", "[SupportSpotsGenerator]") { + const float length = 10; + const float width = 20; + const Polyline polyline{scaled(Vec2f{-length/2.0f, 0.0f}), scaled(Vec2f{length/2.0f, 0.0f})}; + + const Integrals integrals{{polyline}, {width}}; + CHECK(integrals.area == Approx(length * width)); + CHECK(integrals.x_i.x() == Approx(0)); + CHECK(integrals.x_i.y() == Approx(0)); + CHECK(integrals.x_i_squared.x() == Approx(std::pow(length, 3) * width / 12)); + CHECK(integrals.x_i_squared.y() == Approx(length * std::pow(width, 3) / 12)); +} + TEST_CASE("Moment values and ratio check.", "[SupportSpotsGenerator]") { const float width = 40; const float height = 2; From e8064cdb14bce9092829c6965f0d58e7d34f7a43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20=C5=A0ach?= Date: Mon, 20 Nov 2023 10:00:44 +0100 Subject: [PATCH 2/2] Refactor: SupportSpotsGenerator integrals implementation improvements. Split Integrals constructor to accept polygon and polygons. Add operator+(Integrals, Integrals). Make sure extrusion extraction for integral calculation throws in invalid case. --- src/libslic3r/SupportSpotsGenerator.cpp | 86 ++++++++++++------- src/libslic3r/SupportSpotsGenerator.hpp | 8 ++ .../test_support_spots_generator.cpp | 39 +++++---- 3 files changed, 86 insertions(+), 47 deletions(-) diff --git a/src/libslic3r/SupportSpotsGenerator.cpp b/src/libslic3r/SupportSpotsGenerator.cpp index 7ca54793ec..58964dc33f 100644 --- a/src/libslic3r/SupportSpotsGenerator.cpp +++ b/src/libslic3r/SupportSpotsGenerator.cpp @@ -154,23 +154,34 @@ void SliceConnection::print_info(const std::string &tag) const std::cout << "covariance: " << covariance << std::endl; } -Integrals::Integrals (const Polygons& polygons) { +Integrals::Integrals(const Polygon &polygon) +{ + if (polygon.points.size() < 3) { + assert(false && "Polygon is expected to have non-zero area!"); + *this = Integrals{}; + return; + } + Vec2f p0 = unscaled(polygon.first_point()).cast(); + for (size_t i = 2; i < polygon.points.size(); i++) { + Vec2f p1 = unscaled(polygon.points[i - 1]).cast(); + Vec2f p2 = unscaled(polygon.points[i]).cast(); + + float sign = cross2(p1 - p0, p2 - p1) > 0 ? 1.0f : -1.0f; + + auto [area, first_moment_of_area, second_moment_area, + second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2); + + this->area += sign * area; + this->x_i += sign * first_moment_of_area; + this->x_i_squared += sign * second_moment_area; + this->xy += sign * second_moment_of_area_covariance; + } +} + +Integrals::Integrals(const Polygons &polygons) +{ for (const Polygon &polygon : polygons) { - Vec2f p0 = unscaled(polygon.first_point()).cast(); - for (size_t i = 2; i < polygon.points.size(); i++) { - Vec2f p1 = unscaled(polygon.points[i - 1]).cast(); - Vec2f p2 = unscaled(polygon.points[i]).cast(); - - float sign = cross2(p1 - p0, p2 - p1) > 0 ? 1.0f : -1.0f; - - auto [area, first_moment_of_area, second_moment_area, - second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2); - - this->area += sign * area; - this->x_i += sign * first_moment_of_area; - this->x_i_squared += sign * second_moment_area; - this->xy += sign * second_moment_of_area_covariance; - } + *this = *this + Integrals{polygon}; } } @@ -192,15 +203,21 @@ Integrals::Integrals(const Polylines& polylines, const std::vector& width Vec2crd d = scaled(Vec2f{line_a - normal * width/2}); const Polygon ractangle({a, b, c, d}); - Integrals integrals{{ractangle}}; - this->area += integrals.area; - this->x_i += integrals.x_i; - this->x_i_squared += integrals.x_i_squared; - this->xy += integrals.xy; + Integrals integrals{ractangle}; + *this = *this + integrals; } } } +Integrals::Integrals(float area, Vec2f x_i, Vec2f x_i_squared, float xy) + : area(area), x_i(std::move(x_i)), x_i_squared(std::move(x_i_squared)), xy(xy) +{} + +Integrals operator+(const Integrals &a, const Integrals &b) +{ + return Integrals{a.area + b.area, a.x_i + b.x_i, a.x_i_squared + b.x_i_squared, a.xy + b.xy}; +} + SliceConnection estimate_slice_connection(size_t slice_idx, const Layer *layer) { SliceConnection connection; @@ -503,26 +520,33 @@ ObjectPart::ObjectPart( Polylines polylines; std::vector widths; - const auto* path = dynamic_cast(entity); - if (path !=nullptr) { + if ( + const auto* path = dynamic_cast(entity); + path != nullptr + ) { polylines.push_back(path->as_polyline()); widths.push_back(path->width()); - } - - const auto* loop = dynamic_cast(entity); - if (loop !=nullptr) { + } else if ( + const auto* loop = dynamic_cast(entity); + loop != nullptr + ) { for (const ExtrusionPath& path : loop->paths) { polylines.push_back(path.as_polyline()); widths.push_back(path.width()); } - } - - const auto* multi_path = dynamic_cast(entity); - if (multi_path !=nullptr) { + } else if ( + const auto* multi_path = dynamic_cast(entity); + multi_path != nullptr + ) { for (const ExtrusionPath& path : multi_path->paths) { polylines.push_back(path.as_polyline()); widths.push_back(path.width()); } + } else { + throw std::runtime_error( + "Failed to construct object part from extrusions!" + " Unknown extrusion type." + ); } const Integrals integrals{polylines, widths}; diff --git a/src/libslic3r/SupportSpotsGenerator.hpp b/src/libslic3r/SupportSpotsGenerator.hpp index 0f08550b7a..07f264e7da 100644 --- a/src/libslic3r/SupportSpotsGenerator.hpp +++ b/src/libslic3r/SupportSpotsGenerator.hpp @@ -151,10 +151,16 @@ class Integrals{ * @param polygons List of polygons specifing the domain. */ explicit Integrals(const Polygons& polygons); + explicit Integrals(const Polygon& polygon); + /** + * Construct integral x_i int x_i^2 (i=1,2), xy and integral 1 (area) over + * a set of rectangles defined by a "thick" polyline. + */ explicit Integrals(const Polylines& polylines, const std::vector& widths); // TODO refactor and delete the default constructor Integrals() = default; + Integrals(float area, Vec2f x_i, Vec2f x_i_squared, float xy); float area{}; Vec2f x_i{Vec2f::Zero()}; @@ -165,6 +171,8 @@ private: void add(const Integrals& other); }; +Integrals operator+(const Integrals& a, const Integrals& b); + float compute_second_moment( const Integrals& integrals, const Vec2f& axis_direction diff --git a/tests/libslic3r/test_support_spots_generator.cpp b/tests/libslic3r/test_support_spots_generator.cpp index 1b51ccb0ac..7b8dd77bae 100644 --- a/tests/libslic3r/test_support_spots_generator.cpp +++ b/tests/libslic3r/test_support_spots_generator.cpp @@ -5,23 +5,31 @@ using namespace Slic3r; using namespace SupportSpotsGenerator; +namespace Rectangle { +const float width = 10; +const float height = 20; +const Polygon polygon = { + scaled(Vec2f{-width / 2, -height / 2}), + scaled(Vec2f{width / 2, -height / 2}), + scaled(Vec2f{width / 2, height / 2}), + scaled(Vec2f{-width / 2, height / 2}) +}; +} TEST_CASE("Numerical integral over polygon calculation compared with exact solution.", "[SupportSpotsGenerator]") { - const float width = 10; - const float height = 20; - const Polygon polygon = { - scaled(Vec2f{-width / 2, -height / 2}), - scaled(Vec2f{width / 2, -height / 2}), - scaled(Vec2f{width / 2, height / 2}), - scaled(Vec2f{-width / 2, height / 2}) - }; + const Integrals integrals{Rectangle::polygon}; - const Integrals integrals{{polygon}}; - CHECK(integrals.area == Approx(width * height)); + CHECK(integrals.area == Approx(Rectangle::width * Rectangle::height)); CHECK(integrals.x_i.x() == Approx(0)); CHECK(integrals.x_i.y() == Approx(0)); - CHECK(integrals.x_i_squared.x() == Approx(std::pow(width, 3) * height / 12)); - CHECK(integrals.x_i_squared.y() == Approx(width * std::pow(height, 3) / 12)); + CHECK(integrals.x_i_squared.x() == Approx(std::pow(Rectangle::width, 3) * Rectangle::height / 12)); + CHECK(integrals.x_i_squared.y() == Approx(Rectangle::width * std::pow(Rectangle::height, 3) / 12)); +} + +TEST_CASE("Integrals over multiple polygons", "[SupportSpotsGenerator]") { + const Integrals integrals{{Rectangle::polygon, Rectangle::polygon}}; + + CHECK(integrals.area == Approx(2 * Rectangle::width * Rectangle::height)); } TEST_CASE("Numerical integral over line calculation compared with exact solution.", "[SupportSpotsGenerator]") { @@ -50,7 +58,7 @@ TEST_CASE("Moment values and ratio check.", "[SupportSpotsGenerator]") { scaled(Vec2f{0, height}) }; - const Integrals integrals{{polygon}}; + const Integrals integrals{polygon}; const Vec2f x_axis{1, 0}; const float x_axis_moment = compute_second_moment(integrals, x_axis); @@ -68,7 +76,6 @@ TEST_CASE("Moment values and ratio check.", "[SupportSpotsGenerator]") { } TEST_CASE("Moments calculation for rotated axis.", "[SupportSpotsGenerator]") { - Polygon polygon = { scaled(Vec2f{6.362284076172198, 138.9674202217155}), scaled(Vec2f{97.48779843751677, 106.08136606617076}), @@ -82,7 +89,7 @@ TEST_CASE("Moments calculation for rotated axis.", "[SupportSpotsGenerator]") { scaled(Vec2f{77.56229640885199, 189.33057746591336}) }; - Integrals integrals{{polygon}}; + Integrals integrals{polygon}; // Meassured counterclockwise from (1, 0) const float angle = 1.432f; @@ -143,7 +150,7 @@ TEST_CASE_METHOD(ObjectPartFixture, "Constructing ObjectPart using extrusion col std::nullopt }; - Integrals expected{{expected_polygon}}; + Integrals expected{expected_polygon}; CHECK(part.connected_to_bed == true); Vec3f volume_centroid{part.volume_centroid_accumulator / part.volume};