From 6164051d6026308bbb798c56efd8887a46fe83cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20=C5=A0ach?= Date: Tue, 19 Sep 2023 14:55:38 +0200 Subject: [PATCH] Use proper formula for second moment of area. SupportSpotsGenerator originally used a heuristic formula. Current formula is properly derived using known properties of second moment of area. Several tests of this formula are added. --- src/libslic3r/SupportSpotsGenerator.cpp | 168 +++++++++++------- src/libslic3r/SupportSpotsGenerator.hpp | 27 +++ tests/libslic3r/CMakeLists.txt | 1 + .../test_support_spots_generator.cpp | 97 ++++++++++ 4 files changed, 229 insertions(+), 64 deletions(-) create mode 100644 tests/libslic3r/test_support_spots_generator.cpp diff --git a/src/libslic3r/SupportSpotsGenerator.cpp b/src/libslic3r/SupportSpotsGenerator.cpp index fd5bd7d5e8..eb3f1f07a1 100644 --- a/src/libslic3r/SupportSpotsGenerator.cpp +++ b/src/libslic3r/SupportSpotsGenerator.cpp @@ -166,7 +166,7 @@ struct SliceConnection this->second_moment_of_area_covariance_accumulator += other.second_moment_of_area_covariance_accumulator; } - void print_info(const std::string &tag) + void print_info(const std::string &tag) const { Vec3f centroid = centroid_accumulator / area; Vec2f variance = (second_moment_of_area_accumulator / area - centroid.head<2>().cwiseProduct(centroid.head<2>())); @@ -179,6 +179,27 @@ struct SliceConnection } }; +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; + } + } +} + + SliceConnection estimate_slice_connection(size_t slice_idx, const Layer *layer) { SliceConnection connection; @@ -200,22 +221,11 @@ SliceConnection estimate_slice_connection(size_t slice_idx, const Layer *layer) Polygons overlap = intersection(ClipperUtils::clip_clipper_polygons_with_subject_bbox(slice_polys, below_bb), ClipperUtils::clip_clipper_polygons_with_subject_bbox(below_polys, slice_bb)); - for (const Polygon &poly : overlap) { - 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 [area, first_moment_of_area, second_moment_area, - second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2); - connection.area += sign * area; - connection.centroid_accumulator += sign * Vec3f(first_moment_of_area.x(), first_moment_of_area.y(), layer->print_z * area); - connection.second_moment_of_area_accumulator += sign * second_moment_area; - connection.second_moment_of_area_covariance_accumulator += sign * second_moment_of_area_covariance; - } - } + const Integrals integrals{overlap}; + connection.area += integrals.area; + connection.centroid_accumulator += Vec3f(integrals.x_i.x(), integrals.x_i.y(), layer->print_z * integrals.area); + connection.second_moment_of_area_accumulator += integrals.x_i_squared; + connection.second_moment_of_area_covariance_accumulator += integrals.xy; return connection; }; @@ -450,6 +460,48 @@ std::vector check_extrusion_entity_stability(const ExtrusionEntit } } +/** + * Calculates the second moment of area over an arbitrary polygon. + * + * Important note: The calculated moment is for an axis with origin at + * the polygon centroid! + * + * @param integrals Integrals over the polygon area. + * @param axis_direction Direction of the rotation axis going through centroid. + */ +float compute_second_moment( + const Integrals& integrals, + const Vec2f& axis_direction +) { + // Second moment of area for any axis intersecting coordinate system origin + // can be evaluated using the second moments of area calculated for the coordinate + // system axis and the moment product (int xy). + // The equation is derived appling known formulas for the moment of inertia + // to a plannar problem. One can reason about second moment + // of area by by setting density to 1 in the moment of inertia formulas. + const auto area = integrals.area; + const auto I_xx = integrals.x_i_squared.y(); + const auto I_yy = integrals.x_i_squared.x(); + const auto I_xy = -integrals.xy; + + const Vec2f centroid = integrals.x_i / area; + + Matrix2f moment_tensor{}; + moment_tensor << + I_xx, I_xy, + I_xy, I_yy; + + const float moment_at_0_0 = axis_direction.transpose() * moment_tensor * axis_direction; + + // Apply parallel axis theorem to move the moment to centroid + using line_alg::distance_to_infinite_squared; + + const Linef axis_at_0_0 = {{0, 0}, axis_direction.cast()}; + + const double distance = distance_to_infinite_squared(axis_at_0_0, centroid.cast()); + return moment_at_0_0 - area * distance; +} + class ObjectPart { public: @@ -482,43 +534,22 @@ public: this->sticking_second_moment_of_area_covariance_accumulator += sticking_area * position.x() * position.y(); } - float compute_directional_xy_variance(const Vec2f &line_dir, - const Vec3f ¢roid_accumulator, - const Vec2f &second_moment_of_area_accumulator, - const float &second_moment_of_area_covariance_accumulator, - const float &area) const - { - assert(area > 0); - Vec3f centroid = centroid_accumulator / area; - Vec2f variance = (second_moment_of_area_accumulator / area - centroid.head<2>().cwiseProduct(centroid.head<2>())); - float covariance = second_moment_of_area_covariance_accumulator / area - centroid.x() * centroid.y(); - // Var(aX+bY)=a^2*Var(X)+b^2*Var(Y)+2*a*b*Cov(X,Y) - float directional_xy_variance = line_dir.x() * line_dir.x() * variance.x() + line_dir.y() * line_dir.y() * variance.y() + - 2.0f * line_dir.x() * line_dir.y() * covariance; -#ifdef DETAILED_DEBUG_LOGS - BOOST_LOG_TRIVIAL(debug) << "centroid: " << centroid.x() << " " << centroid.y() << " " << centroid.z(); - BOOST_LOG_TRIVIAL(debug) << "variance: " << variance.x() << " " << variance.y(); - BOOST_LOG_TRIVIAL(debug) << "covariance: " << covariance; - BOOST_LOG_TRIVIAL(debug) << "directional_xy_variance: " << directional_xy_variance; -#endif - return directional_xy_variance; - } - float compute_elastic_section_modulus(const Vec2f &line_dir, - const Vec3f &extreme_point, - const Vec3f ¢roid_accumulator, - const Vec2f &second_moment_of_area_accumulator, - const float &second_moment_of_area_covariance_accumulator, - const float &area) const - { - float directional_xy_variance = compute_directional_xy_variance(line_dir, centroid_accumulator, second_moment_of_area_accumulator, - second_moment_of_area_covariance_accumulator, area); - if (directional_xy_variance < EPSILON) { return 0.0f; } - Vec3f centroid = centroid_accumulator / area; + float compute_elastic_section_modulus( + const Vec2f &line_dir, + const Vec3f &extreme_point, + const Integrals& integrals + ) const { + float second_moment_of_area = compute_second_moment(integrals, Vec2f{-line_dir.y(), line_dir.x()}); + + if (second_moment_of_area < EPSILON) { return 0.0f; } + + Vec2f centroid = integrals.x_i / integrals.area; float extreme_fiber_dist = line_alg::distance_to(Linef(centroid.head<2>().cast(), (centroid.head<2>() + Vec2f(line_dir.y(), -line_dir.x())).cast()), extreme_point.head<2>().cast()); - float elastic_section_modulus = area * directional_xy_variance / extreme_fiber_dist; + + float elastic_section_modulus = second_moment_of_area / extreme_fiber_dist; #ifdef DETAILED_DEBUG_LOGS BOOST_LOG_TRIVIAL(debug) << "extreme_fiber_dist: " << extreme_fiber_dist; @@ -534,6 +565,12 @@ public: float layer_z, const Params ¶ms) const { + // Note that exteme point is calculated for the current layer, while it should + // be computed for the first layer. The shape of the first layer however changes a lot, + // during support points additions (for organic supports it is not even clear how) + // and during merging. Using the current layer is heuristics and also small optimization, + // as the AABB tree for it is calculated anyways. This heuristic should usually be + // on the safe side. Vec2f line_dir = (extruded_line.b - extruded_line.a).normalized(); const Vec3f &mass_centroid = this->volume_centroid_accumulator / this->volume; float mass = this->volume * params.filament_density; @@ -548,19 +585,19 @@ public: { if (this->sticking_area < EPSILON) return {1.0f, SupportPointCause::UnstableFloatingPart}; + Integrals integrals; + integrals.area = this->sticking_area; + integrals.x_i = this->sticking_centroid_accumulator.head<2>(); + integrals.x_i_squared = this->sticking_second_moment_of_area_accumulator; + integrals.xy = this->sticking_second_moment_of_area_covariance_accumulator; + Vec3f bed_centroid = this->sticking_centroid_accumulator / this->sticking_area; - float bed_yield_torque = -compute_elastic_section_modulus(line_dir, extreme_point, this->sticking_centroid_accumulator, - this->sticking_second_moment_of_area_accumulator, - this->sticking_second_moment_of_area_covariance_accumulator, - this->sticking_area) * - params.get_bed_adhesion_yield_strength(); + float bed_yield_torque = -compute_elastic_section_modulus(line_dir, extreme_point, integrals) * params.get_bed_adhesion_yield_strength(); Vec2f bed_weight_arm = (mass_centroid.head<2>() - bed_centroid.head<2>()); float bed_weight_arm_len = bed_weight_arm.norm(); - float bed_weight_dir_xy_variance = compute_directional_xy_variance(bed_weight_arm, this->sticking_centroid_accumulator, - this->sticking_second_moment_of_area_accumulator, - this->sticking_second_moment_of_area_covariance_accumulator, - this->sticking_area); + + float bed_weight_dir_xy_variance = compute_second_moment(integrals, {-bed_weight_arm.y(), bed_weight_arm.x()}) / this->sticking_area; float bed_weight_sign = bed_weight_arm_len < 2.0f * sqrt(bed_weight_dir_xy_variance) ? -1.0f : 1.0f; float bed_weight_torque = bed_weight_sign * bed_weight_arm_len * weight; @@ -600,11 +637,14 @@ public: Vec3f conn_centroid = connection.centroid_accumulator / connection.area; if (layer_z - conn_centroid.z() < 3.0f) { return {-1.0f, SupportPointCause::WeakObjectPart}; } - float conn_yield_torque = compute_elastic_section_modulus(line_dir, extreme_point, connection.centroid_accumulator, - connection.second_moment_of_area_accumulator, - connection.second_moment_of_area_covariance_accumulator, - connection.area) * - params.material_yield_strength; + + Integrals integrals; + integrals.area = connection.area; + integrals.x_i = connection.centroid_accumulator.head<2>(); + integrals.x_i_squared = connection.second_moment_of_area_accumulator; + integrals.xy = connection.second_moment_of_area_covariance_accumulator; + + float conn_yield_torque = compute_elastic_section_modulus(line_dir, extreme_point, integrals) * params.material_yield_strength; float conn_weight_arm = (conn_centroid.head<2>() - mass_centroid.head<2>()).norm(); if (layer_z - conn_centroid.z() < 30.0) { diff --git a/src/libslic3r/SupportSpotsGenerator.hpp b/src/libslic3r/SupportSpotsGenerator.hpp index 2e07ef0290..ddb87f0f94 100644 --- a/src/libslic3r/SupportSpotsGenerator.hpp +++ b/src/libslic3r/SupportSpotsGenerator.hpp @@ -147,6 +147,33 @@ struct PartialObject bool connected_to_bed; }; + +/** + * Unsacled values of integrals over a polygonal domain. + */ +class Integrals{ + public: + /** + * Construct integral x_i int x_i^2 (i=1,2), xy and integral 1 (area). + * + * @param polygons List of polygons specifing the domain. + */ + explicit Integrals(const Polygons& polygons); + + // TODO refactor and delete the default constructor + Integrals() = default; + + float area{}; + Vec2f x_i{Vec2f::Zero()}; + Vec2f x_i_squared{Vec2f::Zero()}; + float xy{}; +}; + +float compute_second_moment( + const Integrals& integrals, + const Vec2f& axis_direction +); + using PartialObjects = std::vector; // Both support points and partial objects are sorted from the lowest z to the highest diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt index 90693c8a04..7b0d200c7d 100644 --- a/tests/libslic3r/CMakeLists.txt +++ b/tests/libslic3r/CMakeLists.txt @@ -39,6 +39,7 @@ add_executable(${_TEST_NAME}_tests test_astar.cpp test_anyptr.cpp test_jump_point_search.cpp + test_support_spots_generator.cpp ../data/prusaparts.cpp ../data/prusaparts.hpp ) diff --git a/tests/libslic3r/test_support_spots_generator.cpp b/tests/libslic3r/test_support_spots_generator.cpp new file mode 100644 index 0000000000..3ea8dcd070 --- /dev/null +++ b/tests/libslic3r/test_support_spots_generator.cpp @@ -0,0 +1,97 @@ +#include "libslic3r/Point.hpp" +#include +#include + +using namespace Slic3r; +using namespace SupportSpotsGenerator; + + +TEST_CASE("Numerical integral 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{{polygon}}; + CHECK(integrals.area == Approx(width * 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)); +} + +TEST_CASE("Moment values and ratio check.", "[SupportSpotsGenerator]") { + const float width = 40; + const float height = 2; + + // Moments are calculated at centroid. + // Polygon centroid must not be (0, 0). + const Polygon polygon = { + scaled(Vec2f{0, 0}), + scaled(Vec2f{width, 0}), + scaled(Vec2f{width, height}), + scaled(Vec2f{0, height}) + }; + + const Integrals integrals{{polygon}}; + + const Vec2f x_axis{1, 0}; + const float x_axis_moment = compute_second_moment(integrals, x_axis); + + const Vec2f y_axis{0, 1}; + const float y_axis_moment = compute_second_moment(integrals, y_axis); + + const float moment_ratio = std::pow(width / height, 2); + + // Ensure the object transaltion has no effect. + CHECK(x_axis_moment == Approx(width * std::pow(height, 3) / 12)); + CHECK(y_axis_moment == Approx(std::pow(width, 3) * height / 12)); + // If the object is "wide" the y axis moments should be large compared to x axis moment. + CHECK(y_axis_moment / x_axis_moment == Approx(moment_ratio)); +} + +TEST_CASE("Moments calculation for rotated axis.", "[SupportSpotsGenerator]") { + + Polygon polygon = { + scaled(Vec2f{6.362284076172198, 138.9674202217155}), + scaled(Vec2f{97.48779843751677, 106.08136606617076}), + scaled(Vec2f{135.75221821532384, 66.84428834668765}), + scaled(Vec2f{191.5308049852741, 45.77905628725614}), + scaled(Vec2f{182.7525148049201, 74.01799041087513}), + scaled(Vec2f{296.83210979283473, 196.80022572637228}), + scaled(Vec2f{215.16434429179148, 187.45715418834143}), + scaled(Vec2f{64.64574271229334, 284.293883209721}), + scaled(Vec2f{110.76507036894843, 174.35633141113783}), + scaled(Vec2f{77.56229640885199, 189.33057746591336}) + }; + + Integrals integrals{{polygon}}; + + std::mt19937 generator{std::random_device{}()}; + std::uniform_real_distribution angle_distribution{0, 2*M_PI}; + + // Meassured counterclockwise from (1, 0) + const float angle = angle_distribution(generator); + Vec2f axis{std::cos(angle), std::sin(angle)}; + + float moment_calculated_then_rotated = compute_second_moment( + integrals, + axis + ); + + // We want to rotate the object clockwise by angle to align the axis with (1, 0) + // Method .rotate is counterclockwise for positive angle + polygon.rotate(-angle); + + Integrals integrals_rotated{{polygon}}; + float moment_rotated_polygon = compute_second_moment( + integrals_rotated, + Vec2f{1, 0} + ); + + CHECK(moment_calculated_then_rotated == Approx(moment_rotated_polygon)); +}