amd/orcaslicer
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions 1

ENH: improve bridge direction detection

Browse source 
Thanks to Prusa.
Jira: STUDIO-2453

Change-Id: Iadcc59df44d5abc552f5d558a500fd9bcd66d43f
(cherry picked from commit c19156fd037df4231f3e0cb1e9a899c9b7525372)
This commit is contained in:
Arthur 2023-03-10 20:44:07 +08:00 committed by Lane.Wei
parent dbe1f3f5b1
commit 246ff2653d
6 changed files with 238 additions and 2 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
src/libslic3r/BridgeDetector.cpp
View file

@ -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;

61
src/libslic3r/BridgeDetector.hpp
View file

@ -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 <string>
@ -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<Vec2d, double> 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<double>(), 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<double, Vec2d> directions{};
for (const Line &l : floating_edges) {
Vec2d normal = l.normal().cast<double>().normalized();
double quantized_angle = std::ceil(std::atan2(normal.y(),normal.x()) * 1000.0);
directions.emplace(quantized_angle, normal);
}
std::vector<std::pair<Vec2d, double>> 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<double>();
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<double>::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

2
src/libslic3r/CMakeLists.txt
View file

@ -220,6 +220,8 @@ set(lisbslic3r_sources
PresetBundle.hpp
ProjectTask.cpp
ProjectTask.hpp
PrincipalComponents2D.hpp
PrincipalComponents2D.cpp
AppConfig.cpp
AppConfig.hpp
Print.cpp

10
src/libslic3r/LayerRegion.cpp
View file

@ -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]);
}

140
src/libslic3r/PrincipalComponents2D.cpp Normal file
View file

@ -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<float, Vec2f, Vec2f, float> 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<Vec2f, Vec2f> 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<float>();
for (size_t i = 2; i < poly.points.size(); i++) {
Vec2f p1 = unscaled(poly.points[i - 1]).cast<float>();
Vec2f p2 = unscaled(poly.points[i]).cast<float>();
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<Vec2f, Vec2f> 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};
}
}
}

24
src/libslic3r/PrincipalComponents2D.hpp Normal file
View file

@ -0,0 +1,24 @@
#ifndef slic3r_PrincipalComponents2D_hpp_
#define slic3r_PrincipalComponents2D_hpp_
#include "AABBTreeLines.hpp"
#include "BoundingBox.hpp"
#include "libslic3r.h"
#include <vector>
#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<float, Vec2f, Vec2f, float> 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<Vec2f, Vec2f> compute_principal_components(const Polygons &polys);
}
#endif