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

@@ -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<float>();
+    for (size_t i = 2; i < polygon.points.size(); i++) {
+        Vec2f p1 = unscaled(polygon.points[i - 1]).cast<float>();
+        Vec2f p2 = unscaled(polygon.points[i]).cast<float>();
+
+        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<float>();
+        *this = *this + Integrals{polygon};
-        for (size_t i = 2; i < polygon.points.size(); i++) {
-            Vec2f p1 = unscaled(polygon.points[i - 1]).cast<float>();
-            Vec2f p2 = unscaled(polygon.points[i]).cast<float>();
-
-            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;
-        }
     }
 }
 
@@ -192,15 +203,21 @@ Integrals::Integrals(const Polylines& polylines, const std::vector<float>& width
             Vec2crd d = scaled(Vec2f{line_a - normal * width/2});
 
             const Polygon ractangle({a, b, c, d});
-            Integrals integrals{{ractangle}};
+            Integrals integrals{ractangle};
-            this->area += integrals.area;
+            *this = *this + integrals;
-            this->x_i += integrals.x_i;
-            this->x_i_squared += integrals.x_i_squared;
-            this->xy += integrals.xy;
         }
     }
 }
 
+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<float> widths;
 
-            const auto* path = dynamic_cast<const ExtrusionPath*>(entity);
+            if (
-            if (path !=nullptr) {
+                const auto* path = dynamic_cast<const ExtrusionPath*>(entity);
+                path != nullptr
+            ) {
                 polylines.push_back(path->as_polyline());
                 widths.push_back(path->width());
-            }
+            } else if (
-
+                const auto* loop = dynamic_cast<const ExtrusionLoop*>(entity);
-            const auto* loop = dynamic_cast<const ExtrusionLoop*>(entity);
+                loop != nullptr
-            if (loop !=nullptr) {
+            ) {
                 for (const ExtrusionPath& path : loop->paths) {
                     polylines.push_back(path.as_polyline());
                     widths.push_back(path.width());
                 }
-            }
+            } else if (
-
+                const auto* multi_path = dynamic_cast<const ExtrusionMultiPath*>(entity);
-            const auto* multi_path = dynamic_cast<const ExtrusionMultiPath*>(entity);
+                multi_path != nullptr
-            if (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};

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<float>& 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

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 Integrals integrals{Rectangle::polygon};
-    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(Rectangle::width * Rectangle::height));
-    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.x() == Approx(std::pow(Rectangle::width, 3) * Rectangle::height / 12));
-    CHECK(integrals.x_i_squared.y() == Approx(width * std::pow(height, 3) / 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};