This commit is contained in:
Filip Sykala - NTB T15p 2023-12-04 09:41:14 +01:00
commit 1c0fbe7ff4
3 changed files with 146 additions and 39 deletions

View File

@ -154,26 +154,69 @@ 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>();
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>();
*this = *this + Integrals{polygon};
}
}
float sign = cross2(p1 - p0, p2 - p1) > 0 ? 1.0f : -1.0f;
Integrals::Integrals(const Polylines& polylines, const std::vector<float>& 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<float>();
Vec2f normal{line_direction.y(), -line_direction.x()};
normal.normalize();
auto [area, first_moment_of_area, second_moment_area,
second_moment_of_area_covariance] = compute_moments_of_area_of_triangle(p0, p1, p2);
Vec2f line_a = unscaled(line.a).cast<float>();
Vec2f line_b = unscaled(line.b).cast<float>();
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});
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;
const Polygon ractangle({a, b, c, d});
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)
{
@ -473,18 +516,50 @@ ObjectPart::ObjectPart(
continue;
}
const Polygons polygons{collection->polygons_covered_by_width()};
for (const ExtrusionEntity* entity: collection->flatten()) {
Polylines polylines;
std::vector<float> 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;
if (
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);
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);
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."
);
}
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 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;
}
}
}

View File

@ -151,16 +151,28 @@ 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()};
Vec2f x_i_squared{Vec2f::Zero()};
float xy{};
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

View File

@ -5,23 +5,44 @@
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 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})
};
TEST_CASE("Numerical integral over polygon calculation compared with exact solution.", "[SupportSpotsGenerator]") {
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]") {
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]") {
@ -37,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);
@ -55,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}),
@ -69,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;
@ -130,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};