mirror of
https://git.mirrors.martin98.com/https://github.com/prusa3d/PrusaSlicer.git
synced 2025-07-31 07:32:03 +08:00
Merge branch 'master' into fs_svg_SPE-1517
This commit is contained in:
Filip Sykala - NTB T15p
commit
af520bc787
@ -1,9 +1,12 @@
|
||||
min_slic3r_version = 2.6.2-alpha0
|
||||
1.11.0-alpha5 Added new profiles (additional nozzle diameters) for Prusa MINI Input Shaper (Alpha). Arc fitting changed to I J.
|
||||
1.11.0-alpha4 Updated compatible printer conditions for specific filament profiles.
|
||||
1.11.0-alpha3 Added new print profiles for Prusa MINI Input Shaper (Alpha). Updated MK4 IS profiles.
|
||||
1.11.0-alpha2 Added MK3.9 and Prusa MINI Input Shaper (alpha). Enabled binary g-code, arc fitting and QOI/PNG for MINI and MINI IS.
|
||||
1.11.0-alpha1 Updated ramming parameters. Updated start-gcode for XL Multi-Tool. Updated output filename format.
|
||||
1.11.0-alpha0 Binary g-code, arc fitting, QOI/PNG thumbnails, 90degree XL tower, XL specific filament variants.
|
||||
min_slic3r_version = 2.6.1-rc2
|
||||
1.10.0 Added new profiles for Prusa MINI Input Shaper (Alpha). Enabled XL ramming feature. Prusa XL users may need to re-enable filament profiles in configuration wizard.
|
||||
min_slic3r_version = 2.6.0-beta2
|
||||
1.9.10 Added new print profiles for Prusa MINI Input Shaper (Alpha). Updated MK4 IS profiles.
|
||||
1.9.9 Added Original Prusa MK3.9 and Original Prusa MINI/MINI+ Input Shaper (Alpha). FW version notification (5.0.0 final with input shaper). Updated output filename format. Added additional thumbnail resolution.
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -17,7 +17,9 @@ struct GravityKernel {
|
||||
std::optional<Vec2crd> item_sink;
|
||||
Vec2d active_sink;
|
||||
|
||||
GravityKernel(Vec2crd gravity_center) : sink{gravity_center} {}
|
||||
GravityKernel(Vec2crd gravity_center) :
|
||||
sink{gravity_center}, active_sink{unscaled(gravity_center)} {}
|
||||
|
||||
GravityKernel() = default;
|
||||
|
||||
template<class ArrItem>
|
||||
|
||||
@ -502,6 +502,7 @@ set(SLIC3R_SOURCES
|
||||
Arachne/SkeletalTrapezoidationJoint.hpp
|
||||
Arachne/WallToolPaths.hpp
|
||||
Arachne/WallToolPaths.cpp
|
||||
StaticMap.hpp
|
||||
)
|
||||
|
||||
add_library(libslic3r STATIC ${SLIC3R_SOURCES})
|
||||
|
||||
@ -2641,6 +2641,8 @@ static inline bool validate_smooth_path(const GCode::SmoothPath &smooth_path, bo
|
||||
}
|
||||
#endif //NDEBUG
|
||||
|
||||
static constexpr const double min_gcode_segment_length = 0.002;
|
||||
|
||||
std::string GCodeGenerator::extrude_loop(const ExtrusionLoop &loop_src, const GCode::SmoothPathCache &smooth_path_cache, const std::string_view description, double speed)
|
||||
{
|
||||
// Extrude all loops CCW.
|
||||
@ -2659,7 +2661,7 @@ std::string GCodeGenerator::extrude_loop(const ExtrusionLoop &loop_src, const GC
|
||||
// if polyline was shorter than the clipping distance we'd get a null polyline, so
|
||||
// we discard it in that case.
|
||||
if (m_enable_loop_clipping)
|
||||
clip_end(smooth_path, scale_(EXTRUDER_CONFIG(nozzle_diameter)) * LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER);
|
||||
clip_end(smooth_path, scaled<double>(EXTRUDER_CONFIG(nozzle_diameter)) * LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER, scaled<double>(min_gcode_segment_length));
|
||||
|
||||
if (smooth_path.empty())
|
||||
return {};
|
||||
@ -2704,7 +2706,7 @@ std::string GCodeGenerator::extrude_skirt(
|
||||
// if polyline was shorter than the clipping distance we'd get a null polyline, so
|
||||
// we discard it in that case.
|
||||
if (m_enable_loop_clipping)
|
||||
clip_end(smooth_path, scale_(EXTRUDER_CONFIG(nozzle_diameter)) * LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER);
|
||||
clip_end(smooth_path, scale_(EXTRUDER_CONFIG(nozzle_diameter)) * LOOP_CLIPPING_LENGTH_OVER_NOZZLE_DIAMETER, scaled<double>(min_gcode_segment_length));
|
||||
|
||||
if (smooth_path.empty())
|
||||
return {};
|
||||
@ -3065,12 +3067,14 @@ std::string GCodeGenerator::_extrude(
|
||||
comment = description;
|
||||
comment += description_bridge;
|
||||
}
|
||||
Vec2d prev = this->point_to_gcode_quantized(path.front().point);
|
||||
Vec2d prev_exact = this->point_to_gcode(path.front().point);
|
||||
Vec2d prev = GCodeFormatter::quantize(prev_exact);
|
||||
auto it = path.begin();
|
||||
auto end = path.end();
|
||||
const bool emit_radius = m_config.arc_fitting == ArcFittingType::EmitRadius;
|
||||
for (++ it; it != end; ++ it) {
|
||||
Vec2d p = this->point_to_gcode_quantized(it->point);
|
||||
Vec2d p_exact = this->point_to_gcode(it->point);
|
||||
Vec2d p = GCodeFormatter::quantize(p_exact);
|
||||
assert(p != prev);
|
||||
if (p != prev) {
|
||||
// Center of the radius to be emitted into the G-code: Either by radius or by center offset.
|
||||
@ -3083,11 +3087,12 @@ std::string GCodeGenerator::_extrude(
|
||||
radius = unscaled<double>(it->radius);
|
||||
if (emit_radius) {
|
||||
// Only quantize radius if emitting it directly into G-code. Otherwise use the exact radius for calculating the IJ values.
|
||||
//FIXME rather re-fit the arc to improve accuracy!
|
||||
radius = GCodeFormatter::quantize_xyzf(radius);
|
||||
} else {
|
||||
// Calculate quantized IJ circle center offset.
|
||||
ij = GCodeFormatter::quantize(Vec2d(
|
||||
Geometry::ArcWelder::arc_center(prev.cast<double>(), p.cast<double>(), double(radius), it->ccw())
|
||||
Geometry::ArcWelder::arc_center(prev_exact.cast<double>(), p_exact.cast<double>(), double(radius), it->ccw())
|
||||
- prev));
|
||||
if (ij == Vec2d::Zero())
|
||||
// Don't extrude a degenerated circle.
|
||||
@ -3112,6 +3117,7 @@ std::string GCodeGenerator::_extrude(
|
||||
m_writer.extrude_to_xy_G2G3IJ(p, ij, it->ccw(), dE, comment);
|
||||
}
|
||||
prev = p;
|
||||
prev_exact = p_exact;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,6 +1,9 @@
|
||||
#ifndef slic3r_GCode_LabelObjects_hpp_
|
||||
#define slic3r_GCode_LabelObjects_hpp_
|
||||
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
|
||||
namespace Slic3r {
|
||||
|
||||
enum GCodeFlavor : unsigned char;
|
||||
|
||||
@ -102,7 +102,10 @@ std::optional<Point> sample_path_point_at_distance_from_end(const SmoothPath &pa
|
||||
return {};
|
||||
}
|
||||
|
||||
double clip_end(SmoothPath &path, double distance)
|
||||
// Clip length of a smooth path, for seam hiding.
|
||||
// When clipping the end of a path, don't create segments shorter than min_point_distance_threshold,
|
||||
// rather discard such a degenerate segment.
|
||||
double clip_end(SmoothPath &path, double distance, double min_point_distance_threshold)
|
||||
{
|
||||
while (! path.empty() && distance > 0) {
|
||||
Geometry::ArcWelder::Path &p = path.back().path;
|
||||
@ -111,9 +114,18 @@ double clip_end(SmoothPath &path, double distance)
|
||||
path.pop_back();
|
||||
} else {
|
||||
// Trailing path was trimmed and it is valid.
|
||||
assert(path.back().path.size() > 1);
|
||||
Geometry::ArcWelder::Path &last_path = path.back().path;
|
||||
assert(last_path.size() > 1);
|
||||
assert(distance == 0);
|
||||
// Distance to go is zero.
|
||||
// Remove the last segment if its length is shorter than min_point_distance_threshold.
|
||||
const Geometry::ArcWelder::Segment &prev_segment = last_path[last_path.size() - 2];
|
||||
const Geometry::ArcWelder::Segment &last_segment = last_path.back();
|
||||
if (Geometry::ArcWelder::segment_length<double>(prev_segment, last_segment) < min_point_distance_threshold) {
|
||||
last_path.pop_back();
|
||||
if (last_path.size() < 2)
|
||||
path.pop_back();
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
@ -29,7 +29,9 @@ std::optional<Point> sample_path_point_at_distance_from_start(const SmoothPath &
|
||||
std::optional<Point> sample_path_point_at_distance_from_end(const SmoothPath &path, double distance);
|
||||
|
||||
// Clip end of a smooth path, for seam hiding.
|
||||
double clip_end(SmoothPath &path, double distance);
|
||||
// When clipping the end of a path, don't create segments shorter than min_point_distance_threshold,
|
||||
// rather discard such a degenerate segment.
|
||||
double clip_end(SmoothPath &path, double distance, double min_point_distance_threshold);
|
||||
|
||||
class SmoothPathCache
|
||||
{
|
||||
|
||||
@ -148,36 +148,205 @@ static inline bool get_deviation_sum_squared(const Circle &circle, const Points:
|
||||
return true;
|
||||
}
|
||||
|
||||
double arc_fit_variance(const Point &start_pos, const Point &end_pos, const float radius, bool is_ccw,
|
||||
const Points::const_iterator begin, const Points::const_iterator end)
|
||||
{
|
||||
const Vec2d center = arc_center(start_pos.cast<double>(), end_pos.cast<double>(), double(radius), is_ccw);
|
||||
const double r = std::abs(radius);
|
||||
|
||||
// The circle was calculated from the 1st and last point of the point sequence, thus the fitting of those points does not need to be evaluated.
|
||||
assert(std::abs((begin->cast<double>() - center).norm() - r) < SCALED_EPSILON);
|
||||
assert(std::abs((std::prev(end)->cast<double>() - center).norm() - r) < SCALED_EPSILON);
|
||||
assert(end - begin >= 3);
|
||||
|
||||
double total_deviation = 0;
|
||||
size_t cnt = 0;
|
||||
for (auto it = begin; std::next(it) != end; ++ it) {
|
||||
if (it != begin) {
|
||||
total_deviation += sqr((it->cast<double>() - center).norm() - r);
|
||||
++ cnt;
|
||||
}
|
||||
Point closest_point;
|
||||
if (foot_pt_on_segment(*it, *std::next(it), center.cast<coord_t>(), closest_point)) {
|
||||
total_deviation += sqr((closest_point.cast<double>() - center).cast<double>().norm() - r);
|
||||
++ cnt;
|
||||
}
|
||||
}
|
||||
|
||||
return total_deviation / double(cnt);
|
||||
}
|
||||
|
||||
double arc_fit_max_deviation(const Point &start_pos, const Point &end_pos, const float radius, bool is_ccw,
|
||||
const Points::const_iterator begin, const Points::const_iterator end)
|
||||
{
|
||||
const Vec2d center = arc_center(start_pos.cast<double>(), end_pos.cast<double>(), double(radius), is_ccw);
|
||||
const double r = std::abs(radius);
|
||||
|
||||
// The circle was calculated from the 1st and last point of the point sequence, thus the fitting of those points does not need to be evaluated.
|
||||
assert(std::abs((begin->cast<double>() - center).norm() - r) < SCALED_EPSILON);
|
||||
assert(std::abs((std::prev(end)->cast<double>() - center).norm() - r) < SCALED_EPSILON);
|
||||
assert(end - begin >= 3);
|
||||
|
||||
double max_deviation = 0;
|
||||
double max_signed_deviation = 0;
|
||||
for (auto it = begin; std::next(it) != end; ++ it) {
|
||||
if (it != begin) {
|
||||
double signed_deviation = (it->cast<double>() - center).norm() - r;
|
||||
double deviation = std::abs(signed_deviation);
|
||||
if (deviation > max_deviation) {
|
||||
max_deviation = deviation;
|
||||
max_signed_deviation = signed_deviation;
|
||||
}
|
||||
}
|
||||
Point closest_point;
|
||||
if (foot_pt_on_segment(*it, *std::next(it), center.cast<coord_t>(), closest_point)) {
|
||||
double signed_deviation = (closest_point.cast<double>() - center).cast<double>().norm() - r;
|
||||
double deviation = std::abs(signed_deviation);
|
||||
if (deviation > max_deviation) {
|
||||
max_deviation = deviation;
|
||||
max_signed_deviation = signed_deviation;
|
||||
}
|
||||
}
|
||||
}
|
||||
return max_signed_deviation;
|
||||
}
|
||||
|
||||
static inline int sign(const int64_t i)
|
||||
{
|
||||
return i > 0 ? 1 : i < 0 ? -1 : 0;
|
||||
}
|
||||
|
||||
static std::optional<Circle> try_create_circle(const Points::const_iterator begin, const Points::const_iterator end, const double max_radius, const double tolerance)
|
||||
{
|
||||
std::optional<Circle> out;
|
||||
size_t size = end - begin;
|
||||
if (size == 3) {
|
||||
// Fit the circle throuh the three input points.
|
||||
out = try_create_circle(*begin, *std::next(begin), *std::prev(end), max_radius);
|
||||
if (out && ! circle_approximation_sufficient(*out, begin, end, tolerance))
|
||||
out.reset();
|
||||
} else {
|
||||
#if 0
|
||||
size_t ipivot = size / 2;
|
||||
// Take a center difference of points at the center of the path.
|
||||
//FIXME does it really help? For short arches, the linear interpolation may be
|
||||
Point pivot = (size % 2 == 0) ? (*(begin + ipivot) + *(begin + ipivot - 1)) / 2 :
|
||||
(*(begin + ipivot - 1) + *(begin + ipivot + 1)) / 2;
|
||||
if (std::optional<Circle> circle = try_create_circle(*begin, pivot, *std::prev(end), max_radius);
|
||||
circle && circle_approximation_sufficient(*circle, begin, end, tolerance))
|
||||
return circle;
|
||||
#endif
|
||||
// Find the circle with the least deviation, if one exists.
|
||||
double least_deviation = std::numeric_limits<double>::max();
|
||||
double current_deviation;
|
||||
for (auto it = std::next(begin); std::next(it) != end; ++ it)
|
||||
if (std::optional<Circle> circle = try_create_circle(*begin, *it, *std::prev(end), max_radius);
|
||||
circle && get_deviation_sum_squared(*circle, begin, end, tolerance, current_deviation)) {
|
||||
if (current_deviation < least_deviation) {
|
||||
if (out) {
|
||||
// Fit the center point and the two center points of the two edges with non-linear least squares.
|
||||
std::array<Vec2d, 3> fpts;
|
||||
Vec2d center_point = out->center.cast<double>();
|
||||
Vec2d first_point = begin->cast<double>();
|
||||
Vec2d mid_point = std::next(begin)->cast<double>();
|
||||
Vec2d last_point = std::prev(end)->cast<double>();
|
||||
fpts[0] = 0.5 * (first_point + mid_point);
|
||||
fpts[1] = mid_point;
|
||||
fpts[2] = 0.5 * (mid_point + last_point);
|
||||
const double radius = (first_point - center_point).norm();
|
||||
if (std::abs((fpts[0] - center_point).norm() - radius) < 2. * tolerance &&
|
||||
std::abs((fpts[2] - center_point).norm() - radius) < 2. * tolerance) {
|
||||
if (std::optional<Vec2d> opt_center = ArcWelder::arc_fit_center_gauss_newton_ls(first_point, last_point,
|
||||
center_point, fpts.begin(), fpts.end(), 3);
|
||||
opt_center) {
|
||||
out->center = opt_center->cast<coord_t>();
|
||||
out->radius = (out->radius > 0 ? 1.f : -1.f) * (*opt_center - first_point).norm();
|
||||
}
|
||||
if (! circle_approximation_sufficient(*out, begin, end, tolerance))
|
||||
out.reset();
|
||||
} else
|
||||
out.reset();
|
||||
}
|
||||
} else {
|
||||
std::optional<Circle> circle;
|
||||
{
|
||||
// Try to fit a circle to first, middle and last point.
|
||||
auto mid = begin + (end - begin) / 2;
|
||||
circle = try_create_circle(*begin, *mid, *std::prev(end), max_radius);
|
||||
if (// Use twice the tolerance for fitting the initial circle.
|
||||
// Early exit if such approximation is grossly inaccurate, thus the tolerance could not be achieved.
|
||||
circle && ! circle_approximation_sufficient(*circle, begin, end, tolerance * 2))
|
||||
circle.reset();
|
||||
}
|
||||
if (! circle) {
|
||||
// Find an intersection point of the polyline to be fitted with the bisector of the arc chord.
|
||||
// At such a point the distance of a polyline to an arc wrt. the circle center (or circle radius) will have a largest gradient
|
||||
// of all points on the polyline to be fitted.
|
||||
Vec2i64 first_point = begin->cast<int64_t>();
|
||||
Vec2i64 last_point = std::prev(end)->cast<int64_t>();
|
||||
Vec2i64 c = (first_point.cast<int64_t>() + last_point.cast<int64_t>()) / 2;
|
||||
Vec2i64 v = last_point - first_point;
|
||||
Vec2i64 prev_point = first_point;
|
||||
int prev_side = sign(v.dot(prev_point - c));
|
||||
assert(prev_side != 0);
|
||||
Point point_on_bisector;
|
||||
#ifndef NDEBUG
|
||||
point_on_bisector = { std::numeric_limits<coord_t>::max(), std::numeric_limits<coord_t>::max() };
|
||||
#endif // NDEBUG
|
||||
for (auto it = std::next(begin); it != end; ++ it) {
|
||||
Vec2i64 this_point = it->cast<int64_t>();
|
||||
int64_t d = v.dot(this_point - c);
|
||||
int this_side = sign(d);
|
||||
int sideness = this_side * prev_side;
|
||||
if (sideness < 0) {
|
||||
// Calculate the intersection point.
|
||||
Vec2d vd = v.cast<double>();
|
||||
Vec2d p = c.cast<double>() + vd * double(d) / vd.squaredNorm();
|
||||
point_on_bisector = p.cast<coord_t>();
|
||||
break;
|
||||
}
|
||||
if (sideness == 0) {
|
||||
// this_point is on the bisector.
|
||||
assert(prev_side != 0);
|
||||
assert(this_side == 0);
|
||||
point_on_bisector = this_point.cast<coord_t>();
|
||||
break;
|
||||
}
|
||||
prev_point = this_point;
|
||||
prev_side = this_side;
|
||||
}
|
||||
// point_on_bisector must be set
|
||||
assert(point_on_bisector.x() != std::numeric_limits<coord_t>::max() && point_on_bisector.y() != std::numeric_limits<coord_t>::max());
|
||||
circle = try_create_circle(*begin, point_on_bisector, *std::prev(end), max_radius);
|
||||
if (// Use twice the tolerance for fitting the initial circle.
|
||||
// Early exit if such approximation is grossly inaccurate, thus the tolerance could not be achieved.
|
||||
circle && ! circle_approximation_sufficient(*circle, begin, end, tolerance * 2))
|
||||
circle.reset();
|
||||
}
|
||||
if (circle) {
|
||||
// Fit the arc between the end points by least squares.
|
||||
// Optimize over all points along the path and the centers of the segments.
|
||||
std::vector<Vec2d> fpts;
|
||||
Vec2d first_point = begin->cast<double>();
|
||||
Vec2d last_point = std::prev(end)->cast<double>();
|
||||
Vec2d prev_point = first_point;
|
||||
for (auto it = std::next(begin); it != std::prev(end); ++ it) {
|
||||
Vec2d this_point = it->cast<double>();
|
||||
fpts.emplace_back(0.5 * (prev_point + this_point));
|
||||
fpts.emplace_back(this_point);
|
||||
prev_point = this_point;
|
||||
}
|
||||
fpts.emplace_back(0.5 * (prev_point + last_point));
|
||||
std::optional<Vec2d> opt_center = ArcWelder::arc_fit_center_gauss_newton_ls(first_point, last_point,
|
||||
circle->center.cast<double>(), fpts.begin(), fpts.end(), 5);
|
||||
if (opt_center) {
|
||||
circle->center = opt_center->cast<coord_t>();
|
||||
circle->radius = (circle->radius > 0 ? 1.f : -1.f) * (*opt_center - first_point).norm();
|
||||
if (circle_approximation_sufficient(*circle, begin, end, tolerance)) {
|
||||
out = circle;
|
||||
least_deviation = current_deviation;
|
||||
} else {
|
||||
//FIXME One may consider adjusting the arc to fit the worst offender as a last effort,
|
||||
// however Vojtech is not sure whether it is worth it.
|
||||
}
|
||||
}
|
||||
}
|
||||
/*
|
||||
// From the original arc welder.
|
||||
// Such a loop makes the time complexity of the arc fitting an ugly O(n^3).
|
||||
else {
|
||||
// Find the circle with the least deviation, if one exists.
|
||||
double least_deviation = std::numeric_limits<double>::max();
|
||||
double current_deviation;
|
||||
for (auto it = std::next(begin); std::next(it) != end; ++ it)
|
||||
if (std::optional<Circle> circle = try_create_circle(*begin, *it, *std::prev(end), max_radius);
|
||||
circle && get_deviation_sum_squared(*circle, begin, end, tolerance, current_deviation)) {
|
||||
if (current_deviation < least_deviation) {
|
||||
out = circle;
|
||||
least_deviation = current_deviation;
|
||||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
}
|
||||
return out;
|
||||
}
|
||||
@ -192,11 +361,6 @@ public:
|
||||
Orientation direction { Orientation::Unknown };
|
||||
};
|
||||
|
||||
static inline int sign(const int64_t i)
|
||||
{
|
||||
return i > 0 ? 1 : i < 0 ? -1 : 0;
|
||||
}
|
||||
|
||||
// Return orientation of a polyline with regard to the center.
|
||||
// Successive points are expected to take less than a PI angle step.
|
||||
Orientation arc_orientation(
|
||||
@ -317,21 +481,25 @@ static inline Segments::iterator douglas_peucker_in_place(Segments::iterator beg
|
||||
return douglas_peucker<int64_t>(begin, end, begin, tolerance, [](const Segment &s) { return s.point; });
|
||||
}
|
||||
|
||||
Path fit_path(const Points &src, double tolerance, double fit_circle_percent_tolerance)
|
||||
Path fit_path(const Points &src_in, double tolerance, double fit_circle_percent_tolerance)
|
||||
{
|
||||
assert(tolerance >= 0);
|
||||
assert(fit_circle_percent_tolerance >= 0);
|
||||
double tolerance2 = Slic3r::sqr(tolerance);
|
||||
|
||||
Path out;
|
||||
out.reserve(src.size());
|
||||
if (tolerance <= 0 || src.size() <= 2) {
|
||||
out.reserve(src_in.size());
|
||||
if (tolerance <= 0 || src_in.size() <= 2) {
|
||||
// No simplification, just convert.
|
||||
std::transform(src.begin(), src.end(), std::back_inserter(out), [](const Point &p) -> Segment { return { p }; });
|
||||
} else if (fit_circle_percent_tolerance <= 0) {
|
||||
std::transform(src_in.begin(), src_in.end(), std::back_inserter(out), [](const Point &p) -> Segment { return { p }; });
|
||||
} else if (double tolerance_fine = std::max(0.03 * tolerance, scaled<double>(0.000060));
|
||||
fit_circle_percent_tolerance <= 0 || tolerance_fine > 0.5 * tolerance) {
|
||||
// Convert and simplify to a polyline.
|
||||
std::transform(src.begin(), src.end(), std::back_inserter(out), [](const Point &p) -> Segment { return { p }; });
|
||||
std::transform(src_in.begin(), src_in.end(), std::back_inserter(out), [](const Point &p) -> Segment { return { p }; });
|
||||
out.erase(douglas_peucker_in_place(out.begin(), out.end(), tolerance), out.end());
|
||||
} else {
|
||||
// Simplify the polyline first using a fine threshold.
|
||||
Points src = douglas_peucker(src_in, tolerance_fine);
|
||||
// Perform simplification & fitting.
|
||||
// Index of the start of a last polyline, which has not yet been decimated.
|
||||
int begin_pl_idx = 0;
|
||||
@ -349,13 +517,89 @@ Path fit_path(const Points &src, double tolerance, double fit_circle_percent_tol
|
||||
ArcWelder::default_scaled_max_radius,
|
||||
tolerance, fit_circle_percent_tolerance);
|
||||
this_arc) {
|
||||
// Could extend the arc by one point.
|
||||
assert(this_arc->direction != Orientation::Unknown);
|
||||
arc = this_arc;
|
||||
end = next_end;
|
||||
} else
|
||||
if (end == src.end())
|
||||
// No way to extend the arc.
|
||||
goto fit_end;
|
||||
// Now try to expand the arc by adding points one by one. That should be cheaper than a full arc fit test.
|
||||
while (std::next(end) != src.end()) {
|
||||
assert(end == next_end);
|
||||
{
|
||||
Vec2i64 v1 = arc->start_point.cast<int64_t>() - arc->center.cast<int64_t>();
|
||||
Vec2i64 v2 = arc->end_point.cast<int64_t>() - arc->center.cast<int64_t>();
|
||||
do {
|
||||
if (std::abs((arc->center.cast<double>() - next_end->cast<double>()).norm() - arc->radius) >= tolerance ||
|
||||
inside_arc_wedge_vectors(v1, v2,
|
||||
arc->radius > 0, arc->direction == Orientation::CCW,
|
||||
next_end->cast<int64_t>() - arc->center.cast<int64_t>()))
|
||||
// Cannot extend the current arc with this new point.
|
||||
break;
|
||||
} while (++ next_end != src.end());
|
||||
}
|
||||
if (next_end == end)
|
||||
// No additional point could be added to a current arc.
|
||||
break;
|
||||
// Try to fit a new arc to the extended set of points.
|
||||
// last_tested_failed set to invalid value, no test failed yet.
|
||||
auto last_tested_failed = src.begin();
|
||||
for (;;) {
|
||||
this_arc = try_create_arc(
|
||||
begin, next_end,
|
||||
ArcWelder::default_scaled_max_radius,
|
||||
tolerance, fit_circle_percent_tolerance);
|
||||
if (this_arc) {
|
||||
arc = this_arc;
|
||||
end = next_end;
|
||||
if (last_tested_failed == src.begin()) {
|
||||
// First run of the loop, the arc was extended fully.
|
||||
if (end == src.end())
|
||||
goto fit_end;
|
||||
// Otherwise try to extend the arc with another sample.
|
||||
break;
|
||||
}
|
||||
} else
|
||||
last_tested_failed = next_end;
|
||||
// Take half of the interval up to the failed point.
|
||||
next_end = end + (last_tested_failed - end) / 2;
|
||||
if (next_end == end)
|
||||
// Backed to the last successfull sample.
|
||||
goto fit_end;
|
||||
// Otherwise try to extend the arc up to next_end in another iteration.
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// The last arc was the best we could get.
|
||||
break;
|
||||
}
|
||||
}
|
||||
fit_end:
|
||||
#if 1
|
||||
if (arc) {
|
||||
// Check whether the arc end points are not too close with the risk of quantizing the arc ends to the same point on G-code export.
|
||||
if ((arc->end_point - arc->start_point).cast<double>().squaredNorm() < 2. * sqr(scaled<double>(0.0011))) {
|
||||
// Arc is too short. Skip it, decimate a polyline instead.
|
||||
arc.reset();
|
||||
} else {
|
||||
// Test whether the arc is so flat, that it could be replaced with a straight segment.
|
||||
Line line(arc->start_point, arc->end_point);
|
||||
bool arc_valid = false;
|
||||
for (auto it2 = std::next(begin); it2 != std::prev(end); ++ it2)
|
||||
if (line_alg::distance_to_squared(line, *it2) > tolerance2) {
|
||||
// Polyline could not be fitted by a line segment, thus the arc is considered valid.
|
||||
arc_valid = true;
|
||||
break;
|
||||
}
|
||||
if (! arc_valid)
|
||||
// Arc should be fitted by a line segment. Skip it, decimate a polyline instead.
|
||||
arc.reset();
|
||||
}
|
||||
}
|
||||
#endif
|
||||
if (arc) {
|
||||
// printf("Arc radius: %lf, length: %lf\n", unscaled<double>(arc->radius), arc_length(arc->start_point.cast<double>(), arc->end_point.cast<double>(), arc->radius));
|
||||
// If there is a trailing polyline, decimate it first before saving a new arc.
|
||||
if (out.size() - begin_pl_idx > 2) {
|
||||
// Decimating linear segmens only.
|
||||
@ -363,6 +607,14 @@ Path fit_path(const Points &src, double tolerance, double fit_circle_percent_tol
|
||||
out.erase(douglas_peucker_in_place(out.begin() + begin_pl_idx, out.end(), tolerance), out.end());
|
||||
assert(out.back().linear());
|
||||
}
|
||||
#ifndef NDEBUG
|
||||
// Check for a very short linear segment, that connects two arches. Such segment should not be created.
|
||||
if (out.size() - begin_pl_idx > 1) {
|
||||
const Point& p1 = out[begin_pl_idx].point;
|
||||
const Point& p2 = out.back().point;
|
||||
assert((p2 - p1).cast<double>().squaredNorm() > sqr(scaled<double>(0.0011)));
|
||||
}
|
||||
#endif
|
||||
// Save the index of an end of the new circle segment, which may become the first point of a possible future polyline.
|
||||
begin_pl_idx = int(out.size());
|
||||
// This will be the next point to try to add.
|
||||
@ -413,7 +665,7 @@ Path fit_path(const Points &src, double tolerance, double fit_circle_percent_tol
|
||||
|
||||
#if 0
|
||||
// Verify that all the source points are at tolerance distance from the interpolated path.
|
||||
for (auto it = std::next(src.begin()); it != src.end(); ++ it) {
|
||||
for (auto it = std::next(src_in.begin()); it != src_in.end(); ++ it) {
|
||||
Point start = *std::prev(it);
|
||||
Point end = *it;
|
||||
Vec2d v = (end - start).cast<double>();
|
||||
|
||||
@ -7,7 +7,7 @@
|
||||
|
||||
namespace Slic3r { namespace Geometry { namespace ArcWelder {
|
||||
|
||||
// Calculate center point of an arc given two points and a radius.
|
||||
// Calculate center point (center of a circle) of an arc given two points and a radius.
|
||||
// positive radius: take shorter arc
|
||||
// negative radius: take longer arc
|
||||
// radius must NOT be zero!
|
||||
@ -36,6 +36,37 @@ inline Eigen::Matrix<Float, 2, 1, Eigen::DontAlign> arc_center(
|
||||
return (radius > Float(0)) == is_ccw ? (mid + vp).eval() : (mid - vp).eval();
|
||||
}
|
||||
|
||||
// Calculate middle sample point (point on an arc) of an arc given two points and a radius.
|
||||
// positive radius: take shorter arc
|
||||
// negative radius: take longer arc
|
||||
// radius must NOT be zero!
|
||||
// Taking a sample at the middle of a convex arc (less than PI) is likely much more
|
||||
// useful than taking a sample at the middle of a concave arc (more than PI).
|
||||
template<typename Derived, typename Derived2, typename Float>
|
||||
inline Eigen::Matrix<Float, 2, 1, Eigen::DontAlign> arc_middle_point(
|
||||
const Eigen::MatrixBase<Derived> &start_pos,
|
||||
const Eigen::MatrixBase<Derived2> &end_pos,
|
||||
const Float radius,
|
||||
const bool is_ccw)
|
||||
{
|
||||
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "arc_center(): first parameter is not a 2D vector");
|
||||
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "arc_center(): second parameter is not a 2D vector");
|
||||
static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value, "arc_center(): Both vectors must be of the same type.");
|
||||
static_assert(std::is_same<typename Derived::Scalar, Float>::value, "arc_center(): Radius must be of the same type as the vectors.");
|
||||
assert(radius != 0);
|
||||
using Vector = Eigen::Matrix<Float, 2, 1, Eigen::DontAlign>;
|
||||
auto v = end_pos - start_pos;
|
||||
Float q2 = v.squaredNorm();
|
||||
assert(q2 > 0);
|
||||
Float t2 = sqr(radius) / q2 - Float(.25f);
|
||||
// If the start_pos and end_pos are nearly antipodal, t2 may become slightly negative.
|
||||
// In that case return a centroid of start_point & end_point.
|
||||
Float t = (t2 > 0 ? sqrt(t2) : Float(0)) - radius / sqrt(q2);
|
||||
auto mid = Float(0.5) * (start_pos + end_pos);
|
||||
Vector vp{ -v.y() * t, v.x() * t };
|
||||
return (radius > Float(0)) == is_ccw ? (mid + vp).eval() : (mid - vp).eval();
|
||||
}
|
||||
|
||||
// Calculate angle of an arc given two points and a radius.
|
||||
// Returned angle is in the range <0, 2 PI)
|
||||
// positive radius: take shorter arc
|
||||
@ -87,7 +118,7 @@ inline typename Derived::Scalar arc_length(
|
||||
{
|
||||
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "arc_length(): first parameter is not a 2D vector");
|
||||
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "arc_length(): second parameter is not a 2D vector");
|
||||
static_assert(Derived3::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "arc_length(): third parameter is not a 2D vector");
|
||||
static_assert(Derived3::IsVectorAtCompileTime && int(Derived3::SizeAtCompileTime) == 2, "arc_length(): third parameter is not a 2D vector");
|
||||
static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value &&
|
||||
std::is_same<typename Derived::Scalar, typename Derived3::Scalar>::value, "arc_length(): All third points must be of the same type.");
|
||||
using Float = typename Derived::Scalar;
|
||||
@ -103,6 +134,145 @@ inline typename Derived::Scalar arc_length(
|
||||
return angle * radius;
|
||||
}
|
||||
|
||||
// Be careful! This version has a strong bias towards small circles with small radii
|
||||
// for small angle (nearly straight) arches!
|
||||
// One should rather use arc_fit_center_gauss_newton_ls(), which solves a non-linear least squares problem.
|
||||
//
|
||||
// Calculate center point (center of a circle) of an arc given two fixed points to interpolate
|
||||
// and an additional list of points to fit by least squares.
|
||||
// The circle fitting problem is non-linear, it was linearized by taking difference of squares of radii as a residual.
|
||||
// Therefore the center point is required as a point to linearize at.
|
||||
// Start & end point must be different and the interpolated points must not be collinear with input points.
|
||||
template<typename Derived, typename Derived2, typename Derived3, typename Iterator>
|
||||
inline typename Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign> arc_fit_center_algebraic_ls(
|
||||
const Eigen::MatrixBase<Derived> &start_pos,
|
||||
const Eigen::MatrixBase<Derived2> &end_pos,
|
||||
const Eigen::MatrixBase<Derived3> ¢er_pos,
|
||||
const Iterator it_begin,
|
||||
const Iterator it_end)
|
||||
{
|
||||
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "arc_fit_center_algebraic_ls(): start_pos is not a 2D vector");
|
||||
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "arc_fit_center_algebraic_ls(): end_pos is not a 2D vector");
|
||||
static_assert(Derived3::IsVectorAtCompileTime && int(Derived3::SizeAtCompileTime) == 2, "arc_fit_center_algebraic_ls(): third parameter is not a 2D vector");
|
||||
static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value &&
|
||||
std::is_same<typename Derived::Scalar, typename Derived3::Scalar>::value, "arc_fit_center_algebraic_ls(): All third points must be of the same type.");
|
||||
using Float = typename Derived::Scalar;
|
||||
using Vector = Eigen::Matrix<Float, 2, 1, Eigen::DontAlign>;
|
||||
// Prepare a vector space to transform the fitting into:
|
||||
// center_pos, dir_x, dir_y
|
||||
Vector v = end_pos - start_pos;
|
||||
Vector c = Float(.5) * (start_pos + end_pos);
|
||||
Float lv = v.norm();
|
||||
assert(lv > 0);
|
||||
Vector dir_y = v / lv;
|
||||
Vector dir_x = perp(dir_y);
|
||||
// Center of the new system:
|
||||
// Center X at the projection of center_pos
|
||||
Float offset_x = dir_x.dot(center_pos);
|
||||
// Center is supposed to lie on bisector of the arc end points.
|
||||
// Center Y at the mid point of v.
|
||||
Float offset_y = dir_y.dot(c);
|
||||
assert(std::abs(dir_y.dot(center_pos) - offset_y) < SCALED_EPSILON);
|
||||
assert((dir_x * offset_x + dir_y * offset_y - center_pos).norm() < SCALED_EPSILON);
|
||||
// Solve the least squares fitting in a transformed space.
|
||||
Float a = Float(0.5) * lv;
|
||||
Float b = c.dot(dir_x) - offset_x;
|
||||
Float ab2 = sqr(a) + sqr(b);
|
||||
Float num = Float(0);
|
||||
Float denom = Float(0);
|
||||
const Float w = it_end - it_begin;
|
||||
for (Iterator it = it_begin; it != it_end; ++ it) {
|
||||
Vector p = *it;
|
||||
Float x_i = dir_x.dot(p) - offset_x;
|
||||
Float y_i = dir_y.dot(p) - offset_y;
|
||||
Float x_i2 = sqr(x_i);
|
||||
Float y_i2 = sqr(y_i);
|
||||
num += (x_i - b) * (x_i2 + y_i2 - ab2);
|
||||
denom += b * (b - Float(2) * x_i) + sqr(x_i) + Float(0.25) * w;
|
||||
}
|
||||
assert(denom != 0);
|
||||
Float c_x = Float(0.5) * num / denom;
|
||||
// Transform the center back.
|
||||
Vector out = dir_x * (c_x + offset_x) + dir_y * offset_y;
|
||||
return out;
|
||||
}
|
||||
|
||||
// Calculate center point (center of a circle) of an arc given two fixed points to interpolate
|
||||
// and an additional list of points to fit by non-linear least squares.
|
||||
// The non-linear least squares problem is solved by a Gauss-Newton iterative method.
|
||||
// Start & end point must be different and the interpolated points must not be collinear with input points.
|
||||
// Center position is used to calculate the initial solution of the Gauss-Newton method.
|
||||
//
|
||||
// In case the input points are collinear or close to collinear (such as a small angle arc),
|
||||
// the solution may not converge and an error is indicated.
|
||||
template<typename Derived, typename Derived2, typename Derived3, typename Iterator>
|
||||
inline std::optional<typename Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign>> arc_fit_center_gauss_newton_ls(
|
||||
const Eigen::MatrixBase<Derived> &start_pos,
|
||||
const Eigen::MatrixBase<Derived2> &end_pos,
|
||||
const Eigen::MatrixBase<Derived3> ¢er_pos,
|
||||
const Iterator it_begin,
|
||||
const Iterator it_end,
|
||||
const size_t num_iterations)
|
||||
{
|
||||
static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "arc_fit_center_gauss_newton_ls(): start_pos is not a 2D vector");
|
||||
static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "arc_fit_center_gauss_newton_ls(): end_pos is not a 2D vector");
|
||||
static_assert(Derived3::IsVectorAtCompileTime && int(Derived3::SizeAtCompileTime) == 2, "arc_fit_center_gauss_newton_ls(): third parameter is not a 2D vector");
|
||||
static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value &&
|
||||
std::is_same<typename Derived::Scalar, typename Derived3::Scalar>::value, "arc_fit_center_gauss_newton_ls(): All third points must be of the same type.");
|
||||
using Float = typename Derived::Scalar;
|
||||
using Vector = Eigen::Matrix<Float, 2, 1, Eigen::DontAlign>;
|
||||
// Prepare a vector space to transform the fitting into:
|
||||
// center_pos, dir_x, dir_y
|
||||
Vector v = end_pos - start_pos;
|
||||
Vector c = Float(.5) * (start_pos + end_pos);
|
||||
Float lv = v.norm();
|
||||
assert(lv > 0);
|
||||
Vector dir_y = v / lv;
|
||||
Vector dir_x = perp(dir_y);
|
||||
// Center is supposed to lie on bisector of the arc end points.
|
||||
// Center Y at the mid point of v.
|
||||
Float offset_y = dir_y.dot(c);
|
||||
// Initial value of the parameter to be optimized iteratively.
|
||||
Float c_x = dir_x.dot(center_pos);
|
||||
// Solve the least squares fitting in a transformed space.
|
||||
Float a = Float(0.5) * lv;
|
||||
Float a2 = sqr(a);
|
||||
Float b = c.dot(dir_x);
|
||||
for (size_t iter = 0; iter < num_iterations; ++ iter) {
|
||||
Float num = Float(0);
|
||||
Float denom = Float(0);
|
||||
Float u = b - c_x;
|
||||
// Current estimate of the circle radius.
|
||||
Float r = sqrt(a2 + sqr(u));
|
||||
assert(r > 0);
|
||||
for (Iterator it = it_begin; it != it_end; ++ it) {
|
||||
Vector p = *it;
|
||||
Float x_i = dir_x.dot(p);
|
||||
Float y_i = dir_y.dot(p) - offset_y;
|
||||
Float v = x_i - c_x;
|
||||
Float y_i2 = sqr(y_i);
|
||||
// Distance of i'th sample from the current circle center.
|
||||
Float r_i = sqrt(sqr(v) + y_i2);
|
||||
if (r_i >= EPSILON) {
|
||||
// Square of residual is differentiable at the current c_x and current sample.
|
||||
// Jacobian: diff(residual, c_x)
|
||||
Float j_i = u / r - v / r_i;
|
||||
num += j_i * (r_i - r);
|
||||
denom += sqr(j_i);
|
||||
} else {
|
||||
// Sample point is on current center of the circle,
|
||||
// therefore the gradient is not defined.
|
||||
}
|
||||
}
|
||||
if (denom == 0)
|
||||
// Fitting diverged, the input points are likely nearly collinear with the arch end points.
|
||||
return std::optional<Vector>();
|
||||
c_x -= num / denom;
|
||||
}
|
||||
// Transform the center back.
|
||||
return std::optional<Vector>(dir_x * c_x + dir_y * offset_y);
|
||||
}
|
||||
|
||||
// Test whether a point is inside a wedge of an arc.
|
||||
template<typename Derived, typename Derived2, typename Derived3>
|
||||
inline bool inside_arc_wedge_vectors(
|
||||
@ -194,6 +364,16 @@ size_t arc_discretization_steps(const FloatType radius, const FloatType angle, c
|
||||
// Returned polygon starts with p1, ends with p2 and it is discretized to guarantee the maximum deviation.
|
||||
Points arc_discretize(const Point &p1, const Point &p2, const double radius, const bool ccw, const double deviation);
|
||||
|
||||
// Variance of the arc fit of points <begin, end).
|
||||
// First and last points of <begin, end) are expected to fit the arc exactly.
|
||||
double arc_fit_variance(const Point &start_point, const Point &end_point, const float radius, bool is_ccw,
|
||||
const Points::const_iterator begin, const Points::const_iterator end);
|
||||
|
||||
// Maximum signed deviation of points <begin, end) (L1) from the arc.
|
||||
// First and last points of <begin, end) are expected to fit the arc exactly.
|
||||
double arc_fit_max_deviation(const Point &start_point, const Point &end_point, const float radius, bool is_ccw,
|
||||
const Points::const_iterator begin, const Points::const_iterator end);
|
||||
|
||||
// 1.2m diameter, maximum given by coord_t
|
||||
static constexpr const double default_scaled_max_radius = scaled<double>(600.);
|
||||
// 0.05mm
|
||||
|
||||
@ -147,7 +147,7 @@ Circled circle_ransac(const Vec2ds& input, size_t iterations, double* min_error)
|
||||
}
|
||||
|
||||
template<typename Solver>
|
||||
Circled circle_least_squares_by_solver(const Vec2ds &input, Solver solver)
|
||||
Circled circle_linear_least_squares_by_solver(const Vec2ds &input, Solver solver)
|
||||
{
|
||||
Circled out;
|
||||
if (input.size() < 3) {
|
||||
@ -172,21 +172,21 @@ Circled circle_least_squares_by_solver(const Vec2ds &input, Solver solver)
|
||||
return out;
|
||||
}
|
||||
|
||||
Circled circle_least_squares_svd(const Vec2ds &input)
|
||||
Circled circle_linear_least_squares_svd(const Vec2ds &input)
|
||||
{
|
||||
return circle_least_squares_by_solver(input,
|
||||
return circle_linear_least_squares_by_solver(input,
|
||||
[](const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic /* 3 */> &A, const Eigen::VectorXd &b)
|
||||
{ return A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b).eval(); });
|
||||
}
|
||||
|
||||
Circled circle_least_squares_qr(const Vec2ds &input)
|
||||
Circled circle_linear_least_squares_qr(const Vec2ds &input)
|
||||
{
|
||||
return circle_least_squares_by_solver(input,
|
||||
return circle_linear_least_squares_by_solver(input,
|
||||
[](const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> &A, const Eigen::VectorXd &b)
|
||||
{ return A.colPivHouseholderQr().solve(b).eval(); });
|
||||
}
|
||||
|
||||
Circled circle_least_squares_normal(const Vec2ds &input)
|
||||
Circled circle_linear_least_squares_normal(const Vec2ds &input)
|
||||
{
|
||||
Circled out;
|
||||
if (input.size() < 3) {
|
||||
|
||||
@ -141,12 +141,17 @@ Circled circle_taubin_newton(const Vec2ds& input, size_t cycles = 20);
|
||||
// Find circle using RANSAC randomized algorithm.
|
||||
Circled circle_ransac(const Vec2ds& input, size_t iterations = 20, double* min_error = nullptr);
|
||||
|
||||
// Least squares fitting with SVD. Most accurate, but slowest.
|
||||
Circled circle_least_squares_svd(const Vec2ds &input);
|
||||
// Least squares fitting with QR decomposition. Medium accuracy, medium speed.
|
||||
Circled circle_least_squares_qr(const Vec2ds &input);
|
||||
// Least squares fitting solving normal equations. Low accuracy, high speed.
|
||||
Circled circle_least_squares_normal(const Vec2ds &input);
|
||||
// Linear Least squares fitting.
|
||||
// Be careful! The linear least squares fitting is strongly biased towards small circles,
|
||||
// thus the method is only recommended for circles or arches with large arc angle.
|
||||
// Also it is strongly recommended to center the input at an expected circle (or arc) center
|
||||
// to minimize the small circle bias!
|
||||
// Linear Least squares fitting with SVD. Most accurate, but slowest.
|
||||
Circled circle_linear_least_squares_svd(const Vec2ds &input);
|
||||
// Linear Least squares fitting with QR decomposition. Medium accuracy, medium speed.
|
||||
Circled circle_linear_least_squares_qr(const Vec2ds &input);
|
||||
// Linear Least squares fitting solving normal equations. Low accuracy, high speed.
|
||||
Circled circle_linear_least_squares_normal(const Vec2ds &input);
|
||||
|
||||
// Randomized algorithm by Emo Welzl, working with squared radii for efficiency. The returned circle radius is inflated by epsilon.
|
||||
template<typename Vector, typename Points>
|
||||
|
||||
328
src/libslic3r/StaticMap.hpp
Normal file
328
src/libslic3r/StaticMap.hpp
Normal file
@ -0,0 +1,328 @@
|
||||
#ifndef PRUSASLICER_STATICMAP_HPP
|
||||
#define PRUSASLICER_STATICMAP_HPP
|
||||
|
||||
#include <optional>
|
||||
#include <array>
|
||||
#include <string_view>
|
||||
|
||||
namespace Slic3r {
|
||||
|
||||
// This module provides std::map and std::set like structures with fixed number
|
||||
// of elements and usable at compile or in time constexpr contexts without
|
||||
// any memory allocations.
|
||||
|
||||
// C++20 emulation utilities to get the missing constexpr functionality in C++17
|
||||
namespace static_set_detail {
|
||||
|
||||
// Simple bubble sort but constexpr
|
||||
template<class T, size_t N, class Cmp = std::less<T>>
|
||||
constexpr void sort_array(std::array<T, N> &arr, Cmp cmp = {})
|
||||
{
|
||||
// A bubble sort will do the job, C++20 will have constexpr std::sort
|
||||
for (size_t i = 0; i < N - 1; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < N - i - 1; ++j)
|
||||
{
|
||||
if (!cmp(arr[j], arr[j + 1])) {
|
||||
T temp = arr[j];
|
||||
arr[j] = arr[j + 1];
|
||||
arr[j + 1] = temp;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Simple emulation of lower_bound with constexpr
|
||||
template<class It, class V, class Cmp = std::less<V>>
|
||||
constexpr auto array_lower_bound(It from, It to, const V &val, Cmp cmp)
|
||||
{
|
||||
auto N = std::distance(from, to);
|
||||
std::size_t middle = N / 2;
|
||||
|
||||
if (N == 0) {
|
||||
return from; // Key not found, return the beginning of the array
|
||||
} else if (cmp(val, *(from + middle))) {
|
||||
return array_lower_bound(from, from + middle, val, cmp);
|
||||
} else if (cmp(*(from + middle), val)) {
|
||||
return array_lower_bound(from + middle + 1, to, val, cmp);
|
||||
} else {
|
||||
return from + middle; // Key found, return an iterator to it
|
||||
}
|
||||
|
||||
return to;
|
||||
}
|
||||
|
||||
template<class T, size_t N, class Cmp = std::less<T>>
|
||||
constexpr auto array_lower_bound(const std::array<T, N> &arr,
|
||||
const T &val,
|
||||
Cmp cmp = {})
|
||||
{
|
||||
return array_lower_bound(arr.begin(), arr.end(), val, cmp);
|
||||
}
|
||||
|
||||
template<class T, std::size_t N, std::size_t... I>
|
||||
constexpr std::array<std::remove_cv_t<T>, N>
|
||||
to_array_impl(T (&a)[N], std::index_sequence<I...>)
|
||||
{
|
||||
return {{a[I]...}};
|
||||
}
|
||||
|
||||
template<class T, std::size_t N>
|
||||
constexpr std::array<std::remove_cv_t<T>, N> to_array(T (&a)[N])
|
||||
{
|
||||
return to_array_impl(a, std::make_index_sequence<N>{});
|
||||
}
|
||||
|
||||
// Emulating constexpr std::pair
|
||||
template<class K, class V>
|
||||
struct StaticMapElement {
|
||||
using KeyType = K;
|
||||
|
||||
constexpr StaticMapElement(const K &k, const V &v): first{k}, second{v} {}
|
||||
|
||||
K first;
|
||||
V second;
|
||||
};
|
||||
|
||||
} // namespace static_set_detail
|
||||
} // namespace Slic3r
|
||||
|
||||
namespace Slic3r {
|
||||
|
||||
// std::set like set structure
|
||||
template<class T, size_t N, class Cmp = std::less<T>>
|
||||
class StaticSet {
|
||||
std::array<T, N> m_vals; // building on top of std::array
|
||||
Cmp m_cmp;
|
||||
|
||||
public:
|
||||
using value_type = T;
|
||||
|
||||
constexpr StaticSet(const std::array<T, N> &arr, Cmp cmp = {})
|
||||
: m_vals{arr}, m_cmp{cmp}
|
||||
{
|
||||
// TODO: C++20 can use std::sort(vals.begin(), vals.end())
|
||||
static_set_detail::sort_array(m_vals, m_cmp);
|
||||
}
|
||||
|
||||
template<class...Ts>
|
||||
constexpr StaticSet(Ts &&...args): m_vals{std::forward<Ts>(args)...}
|
||||
{
|
||||
static_set_detail::sort_array(m_vals, m_cmp);
|
||||
}
|
||||
|
||||
template<class...Ts>
|
||||
constexpr StaticSet(Cmp cmp, Ts &&...args)
|
||||
: m_vals{std::forward<Ts>(args)...}, m_cmp{cmp}
|
||||
{
|
||||
static_set_detail::sort_array(m_vals, m_cmp);
|
||||
}
|
||||
|
||||
constexpr auto find(const T &val) const
|
||||
{
|
||||
// TODO: C++20 can use std::lower_bound
|
||||
auto it = static_set_detail::array_lower_bound(m_vals, val, m_cmp);
|
||||
if (it != m_vals.end() && ! m_cmp(*it, val) && !m_cmp(val, *it) )
|
||||
return it;
|
||||
|
||||
return m_vals.cend();
|
||||
}
|
||||
|
||||
constexpr bool empty() const { return m_vals.empty(); }
|
||||
constexpr size_t size() const { return m_vals.size(); }
|
||||
|
||||
// Can be iterated over
|
||||
constexpr auto begin() const { return m_vals.begin(); }
|
||||
constexpr auto end() const { return m_vals.end(); }
|
||||
};
|
||||
|
||||
// These are "deduction guides", a C++17 feature.
|
||||
// Reason is to be able to deduce template arguments from constructor arguments
|
||||
// e.g.: StaticSet{1, 2, 3} is deduced as StaticSet<int, 3, std::less<int>>, no
|
||||
// need to state the template types explicitly.
|
||||
template<class T, class...Vals>
|
||||
StaticSet(T, Vals...) ->
|
||||
StaticSet<std::enable_if_t<(std::is_same_v<T, Vals> && ...), T>,
|
||||
1 + sizeof...(Vals)>;
|
||||
|
||||
// Same as above, only with the first argument being a comparison functor
|
||||
template<class Cmp, class T, class...Vals>
|
||||
StaticSet(Cmp, T, Vals...) ->
|
||||
StaticSet<std::enable_if_t<(std::is_same_v<T, Vals> && ...), T>,
|
||||
1 + sizeof...(Vals),
|
||||
std::enable_if_t<std::is_invocable_r_v<bool, Cmp, T, T>, Cmp>>;
|
||||
|
||||
// Specialization for the empty set case.
|
||||
template<class T>
|
||||
class StaticSet<T, size_t{0}> {
|
||||
public:
|
||||
constexpr StaticSet() = default;
|
||||
constexpr auto find(const T &val) const { return nullptr; }
|
||||
constexpr bool empty() const { return true; }
|
||||
constexpr size_t size() const { return 0; }
|
||||
constexpr auto begin() const { return nullptr; }
|
||||
constexpr auto end() const { return nullptr; }
|
||||
};
|
||||
|
||||
// Constructor with no arguments need to be deduced as the specialization for
|
||||
// empty sets (see above)
|
||||
StaticSet() -> StaticSet<int, 0>;
|
||||
|
||||
|
||||
|
||||
// StaticMap definition:
|
||||
|
||||
template<class K, class V>
|
||||
using SMapEl = static_set_detail::StaticMapElement<K, V>;
|
||||
|
||||
template<class K, class V>
|
||||
struct DefaultCmp {
|
||||
constexpr bool operator() (const SMapEl<K, V> &el1, const SMapEl<K, V> &el2) const
|
||||
{
|
||||
return std::less<K>{}(el1.first, el2.first);
|
||||
}
|
||||
};
|
||||
|
||||
// Overriding the default comparison for C style strings, as std::less<const char*>
|
||||
// doesn't do the lexicographic comparisons, only the pointer values would be
|
||||
// compared. Fortunately we can wrap the C style strings with string_views and
|
||||
// do the comparison with those.
|
||||
template<class V>
|
||||
struct DefaultCmp<const char *, V> {
|
||||
constexpr bool operator() (const SMapEl<const char *, V> &el1,
|
||||
const SMapEl<const char *, V> &el2) const
|
||||
{
|
||||
return std::string_view{el1.first} < std::string_view{el2.first};
|
||||
}
|
||||
};
|
||||
|
||||
template<class K, class V, size_t N, class Cmp = DefaultCmp<K, V>>
|
||||
class StaticMap {
|
||||
std::array<SMapEl<K, V>, N> m_vals;
|
||||
Cmp m_cmp;
|
||||
|
||||
public:
|
||||
using value_type = SMapEl<K, V>;
|
||||
|
||||
constexpr StaticMap(const std::array<SMapEl<K, V>, N> &arr, Cmp cmp = {})
|
||||
: m_vals{arr}, m_cmp{cmp}
|
||||
{
|
||||
static_set_detail::sort_array(m_vals, cmp);
|
||||
}
|
||||
|
||||
constexpr auto find(const K &key) const
|
||||
{
|
||||
auto ret = m_vals.end();
|
||||
|
||||
SMapEl<K, V> vkey{key, V{}};
|
||||
|
||||
auto it = static_set_detail::array_lower_bound(
|
||||
std::begin(m_vals), std::end(m_vals), vkey, m_cmp
|
||||
);
|
||||
|
||||
if (it != std::end(m_vals) && ! m_cmp(*it, vkey) && !m_cmp(vkey, *it))
|
||||
ret = it;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
constexpr const V& at(const K& key) const
|
||||
{
|
||||
if (auto it = find(key); it != end())
|
||||
return it->second;
|
||||
|
||||
throw std::out_of_range{"No such element"};
|
||||
}
|
||||
|
||||
constexpr bool empty() const { return m_vals.empty(); }
|
||||
constexpr size_t size() const { return m_vals.size(); }
|
||||
|
||||
constexpr auto begin() const { return m_vals.begin(); }
|
||||
constexpr auto end() const { return m_vals.end(); }
|
||||
};
|
||||
|
||||
template<class K, class V>
|
||||
class StaticMap<K, V, size_t{0}> {
|
||||
public:
|
||||
constexpr StaticMap() = default;
|
||||
constexpr auto find(const K &key) const { return nullptr; }
|
||||
constexpr bool empty() const { return true; }
|
||||
constexpr size_t size() const { return 0; }
|
||||
[[noreturn]] constexpr const V& at(const K &) const { throw std::out_of_range{"Map is empty"}; }
|
||||
constexpr auto begin() const { return nullptr; }
|
||||
constexpr auto end() const { return nullptr; }
|
||||
};
|
||||
|
||||
// Deducing template arguments from the StaticMap constructors is not easy,
|
||||
// so there is a helper "make" function to be used instead:
|
||||
// e.g.: auto map = make_staticmap<const char*, int>({ {"one", 1}, {"two", 2}})
|
||||
// will work, and only the key and value type needs to be specified. No need
|
||||
// to state the number of elements, that is deduced automatically.
|
||||
template<class K, class V, size_t N>
|
||||
constexpr auto make_staticmap(const SMapEl<K, V> (&arr) [N])
|
||||
{
|
||||
return StaticMap<K, V, N>{static_set_detail ::to_array(arr), DefaultCmp<K, V>{}};
|
||||
}
|
||||
|
||||
template<class K, class V, size_t N, class Cmp>
|
||||
constexpr auto make_staticmap(const SMapEl<K, V> (&arr) [N], Cmp cmp)
|
||||
{
|
||||
return StaticMap<K, V, N, Cmp>{static_set_detail ::to_array(arr), cmp};
|
||||
}
|
||||
|
||||
// Override for empty maps
|
||||
template<class K, class V, class Cmp = DefaultCmp<K, V>>
|
||||
constexpr auto make_staticmap()
|
||||
{
|
||||
return StaticMap<K, V, 0, Cmp>{};
|
||||
}
|
||||
|
||||
// Override which uses a c++ array as the initializer
|
||||
template<class K, class V, size_t N, class Cmp = DefaultCmp<K, V>>
|
||||
constexpr auto make_staticmap(const std::array<SMapEl<K, V>, N> &arr, Cmp cmp = {})
|
||||
{
|
||||
return StaticMap<K, V, N, Cmp>{arr, cmp};
|
||||
}
|
||||
|
||||
// Helper function to get a specific element from a set, returning a std::optional
|
||||
// which is more convinient than working with iterators
|
||||
template<class V, size_t N, class Cmp, class T>
|
||||
constexpr std::enable_if_t<std::is_convertible_v<T, V>, std::optional<V>>
|
||||
query(const StaticSet<V, N, Cmp> &sset, const T &val)
|
||||
{
|
||||
std::optional<V> ret;
|
||||
if (auto it = sset.find(val); it != sset.end())
|
||||
ret = *it;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
template<class K, class V, size_t N, class Cmp, class KeyT>
|
||||
constexpr std::enable_if_t<std::is_convertible_v<K, KeyT>, std::optional<V>>
|
||||
query(const StaticMap<K, V, N, Cmp> &sset, const KeyT &val)
|
||||
{
|
||||
std::optional<V> ret;
|
||||
|
||||
if (auto it = sset.find(val); it != sset.end())
|
||||
ret = it->second;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
template<class V, size_t N, class Cmp, class T>
|
||||
constexpr std::enable_if_t<std::is_convertible_v<T, V>, bool>
|
||||
contains(const StaticSet<V, N, Cmp> &sset, const T &val)
|
||||
{
|
||||
return sset.find(val) != sset.end();
|
||||
}
|
||||
|
||||
template<class K, class V, size_t N, class Cmp, class KeyT>
|
||||
constexpr std::enable_if_t<std::is_convertible_v<K, KeyT>, bool>
|
||||
contains(const StaticMap<K, V, N, Cmp> &smap, const KeyT &key)
|
||||
{
|
||||
return smap.find(key) != smap.end();
|
||||
}
|
||||
|
||||
} // namespace Slic3r
|
||||
|
||||
#endif // STATICMAP_HPP
|
||||
@ -53,8 +53,12 @@
|
||||
#ifdef DEBUG_FILES
|
||||
#include <boost/nowide/cstdio.hpp>
|
||||
#include "libslic3r/Color.hpp"
|
||||
constexpr bool debug_files = true;
|
||||
#else
|
||||
constexpr bool debug_files = false;
|
||||
#endif
|
||||
|
||||
|
||||
namespace Slic3r::SupportSpotsGenerator {
|
||||
|
||||
ExtrusionLine::ExtrusionLine() : a(Vec2f::Zero()), b(Vec2f::Zero()), len(0.0), origin_entity(nullptr) {}
|
||||
@ -72,10 +76,6 @@ bool ExtrusionLine::is_external_perimeter() const
|
||||
return origin_entity->role().is_external_perimeter();
|
||||
}
|
||||
|
||||
auto get_a(ExtrusionLine &&l) { return l.a; }
|
||||
auto get_b(ExtrusionLine &&l) { return l.b; }
|
||||
|
||||
|
||||
using LD = AABBTreeLines::LinesDistancer<ExtrusionLine>;
|
||||
|
||||
struct SupportGridFilter
|
||||
@ -761,6 +761,284 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
// Function that is used when new support point is generated. It will update the ObjectPart stability, weakest conneciton info,
|
||||
// and the support presence grid and add the point to the issues.
|
||||
void reckon_new_support_point(ObjectPart &part,
|
||||
SliceConnection &weakest_conn,
|
||||
SupportPoints &supp_points,
|
||||
SupportGridFilter &supports_presence_grid,
|
||||
const SupportPoint& support_point,
|
||||
bool is_global = false)
|
||||
{
|
||||
// if position is taken and point is for global stability (force > 0) or we are too close to the bed, do not add
|
||||
// This allows local support points (e.g. bridging) to be generated densely
|
||||
if ((supports_presence_grid.position_taken(support_point.position) && is_global)) {
|
||||
return;
|
||||
}
|
||||
|
||||
float area = support_point.spot_radius * support_point.spot_radius * float(PI);
|
||||
// add the stability effect of the point only if the spot is not taken, so that the densely created local support points do
|
||||
// not add unrealistic amount of stability to the object (due to overlaping of local support points)
|
||||
if (!(supports_presence_grid.position_taken(support_point.position))) {
|
||||
part.add_support_point(support_point.position, area);
|
||||
}
|
||||
|
||||
supp_points.push_back(support_point);
|
||||
supports_presence_grid.take_position(support_point.position);
|
||||
|
||||
// The support point also increases the stability of the weakest connection of the object, which should be reflected
|
||||
if (weakest_conn.area > EPSILON) { // Do not add it to the weakest connection if it is not valid - does not exist
|
||||
weakest_conn.area += area;
|
||||
weakest_conn.centroid_accumulator += support_point.position * area;
|
||||
weakest_conn.second_moment_of_area_accumulator += area *
|
||||
support_point.position.head<2>().cwiseProduct(support_point.position.head<2>());
|
||||
weakest_conn.second_moment_of_area_covariance_accumulator += area * support_point.position.x() * support_point.position.y();
|
||||
}
|
||||
}
|
||||
|
||||
struct LocalSupports {
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> unstable_lines_per_slice;
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> ext_perim_lines_per_slice;
|
||||
};
|
||||
|
||||
struct EnitityToCheck
|
||||
{
|
||||
const ExtrusionEntity *e;
|
||||
const LayerRegion *region;
|
||||
size_t slice_idx;
|
||||
};
|
||||
|
||||
// TODO DRY: Very similar to gather extrusions.
|
||||
std::vector<EnitityToCheck> gather_entities_to_check(const Layer* layer) {
|
||||
auto get_flat_entities = [](const ExtrusionEntity *e) {
|
||||
std::vector<const ExtrusionEntity *> entities;
|
||||
std::vector<const ExtrusionEntity *> queue{e};
|
||||
while (!queue.empty()) {
|
||||
const ExtrusionEntity *next = queue.back();
|
||||
queue.pop_back();
|
||||
if (next->is_collection()) {
|
||||
for (const ExtrusionEntity *e : static_cast<const ExtrusionEntityCollection *>(next)->entities) {
|
||||
queue.push_back(e);
|
||||
}
|
||||
} else {
|
||||
entities.push_back(next);
|
||||
}
|
||||
}
|
||||
return entities;
|
||||
};
|
||||
|
||||
std::vector<EnitityToCheck> entities_to_check;
|
||||
for (size_t slice_idx = 0; slice_idx < layer->lslices_ex.size(); ++slice_idx) {
|
||||
const LayerSlice &slice = layer->lslices_ex.at(slice_idx);
|
||||
for (const auto &island : slice.islands) {
|
||||
for (const LayerExtrusionRange &fill_range : island.fills) {
|
||||
const LayerRegion *fill_region = layer->get_region(fill_range.region());
|
||||
for (size_t fill_idx : fill_range) {
|
||||
for (const ExtrusionEntity *e : get_flat_entities(fill_region->fills().entities[fill_idx])) {
|
||||
if (e->role() == ExtrusionRole::BridgeInfill) {
|
||||
entities_to_check.push_back({e, fill_region, slice_idx});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const LayerRegion *perimeter_region = layer->get_region(island.perimeters.region());
|
||||
for (size_t perimeter_idx : island.perimeters) {
|
||||
for (const ExtrusionEntity *e : get_flat_entities(perimeter_region->perimeters().entities[perimeter_idx])) {
|
||||
entities_to_check.push_back({e, perimeter_region, slice_idx});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return entities_to_check;
|
||||
}
|
||||
|
||||
LocalSupports compute_local_supports(
|
||||
const std::vector<EnitityToCheck>& entities_to_check,
|
||||
const std::optional<Linesf>& previous_layer_boundary,
|
||||
const LD& prev_layer_ext_perim_lines,
|
||||
size_t slices_count,
|
||||
const Params& params
|
||||
) {
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> unstable_lines_per_slice(slices_count);
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> ext_perim_lines_per_slice(slices_count);
|
||||
|
||||
AABBTreeLines::LinesDistancer<Linef> prev_layer_boundary_distancer =
|
||||
(previous_layer_boundary ? AABBTreeLines::LinesDistancer<Linef>{*previous_layer_boundary} : AABBTreeLines::LinesDistancer<Linef>{});
|
||||
|
||||
if constexpr (debug_files) {
|
||||
for (const auto &e_to_check : entities_to_check) {
|
||||
for (const auto &line : check_extrusion_entity_stability(e_to_check.e, e_to_check.region, prev_layer_ext_perim_lines,
|
||||
prev_layer_boundary_distancer, params)) {
|
||||
if (line.support_point_generated.has_value()) {
|
||||
unstable_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
if (line.is_external_perimeter()) {
|
||||
ext_perim_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
tbb::parallel_for(tbb::blocked_range<size_t>(0, entities_to_check.size()),
|
||||
[&entities_to_check, &prev_layer_ext_perim_lines, &prev_layer_boundary_distancer, &unstable_lines_per_slice,
|
||||
&ext_perim_lines_per_slice, ¶ms](tbb::blocked_range<size_t> r) {
|
||||
for (size_t entity_idx = r.begin(); entity_idx < r.end(); ++entity_idx) {
|
||||
const auto &e_to_check = entities_to_check[entity_idx];
|
||||
for (const auto &line :
|
||||
check_extrusion_entity_stability(e_to_check.e, e_to_check.region, prev_layer_ext_perim_lines,
|
||||
prev_layer_boundary_distancer, params)) {
|
||||
if (line.support_point_generated.has_value()) {
|
||||
unstable_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
if (line.is_external_perimeter()) {
|
||||
ext_perim_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
}
|
||||
}
|
||||
});
|
||||
}
|
||||
return {unstable_lines_per_slice, ext_perim_lines_per_slice};
|
||||
}
|
||||
|
||||
struct SliceMappings
|
||||
{
|
||||
std::unordered_map<size_t, size_t> index_to_object_part_mapping;
|
||||
std::unordered_map<size_t, SliceConnection> index_to_weakest_connection;
|
||||
};
|
||||
|
||||
std::optional<PartialObject> to_partial_object(const ObjectPart& part) {
|
||||
if (part.volume > EPSILON) {
|
||||
return PartialObject{part.volume_centroid_accumulator / part.volume, part.volume,
|
||||
part.connected_to_bed};
|
||||
}
|
||||
return {};
|
||||
}
|
||||
|
||||
SliceMappings update_active_object_parts(const Layer *layer,
|
||||
const Params ¶ms,
|
||||
const std::vector<SliceConnection> &precomputed_slice_connections,
|
||||
const SliceMappings &previous_slice_mappings,
|
||||
ActiveObjectParts &active_object_parts,
|
||||
PartialObjects &partial_objects)
|
||||
{
|
||||
SliceMappings new_slice_mappings;
|
||||
|
||||
for (size_t slice_idx = 0; slice_idx < layer->lslices_ex.size(); ++slice_idx) {
|
||||
const LayerSlice &slice = layer->lslices_ex.at(slice_idx);
|
||||
const std::vector<const ExtrusionEntityCollection*> extrusion_collections{gather_extrusions(slice, layer)};
|
||||
const bool connected_to_bed = int(layer->id()) == params.raft_layers_count;
|
||||
|
||||
const std::optional<Polygons> brim{
|
||||
has_brim(layer, params) ?
|
||||
std::optional{get_brim(layer->lslices[slice_idx], params.brim_type, params.brim_width)} :
|
||||
std::nullopt
|
||||
};
|
||||
ObjectPart new_part{
|
||||
extrusion_collections,
|
||||
connected_to_bed,
|
||||
layer->print_z,
|
||||
layer->height,
|
||||
brim
|
||||
};
|
||||
|
||||
const SliceConnection &connection_to_below = precomputed_slice_connections[slice_idx];
|
||||
|
||||
#ifdef DETAILED_DEBUG_LOGS
|
||||
std::cout << "SLICE IDX: " << slice_idx << std::endl;
|
||||
for (const auto &link : slice.overlaps_below) {
|
||||
std::cout << "connected to slice below: " << link.slice_idx << " by area : " << link.area << std::endl;
|
||||
}
|
||||
connection_to_below.print_info("CONNECTION TO BELOW");
|
||||
#endif
|
||||
|
||||
if (connection_to_below.area < EPSILON) { // new object part emerging
|
||||
size_t part_id = active_object_parts.insert(new_part);
|
||||
new_slice_mappings.index_to_object_part_mapping.emplace(slice_idx, part_id);
|
||||
new_slice_mappings.index_to_weakest_connection.emplace(slice_idx, connection_to_below);
|
||||
} else {
|
||||
size_t final_part_id{};
|
||||
SliceConnection transfered_weakest_connection{};
|
||||
// MERGE parts
|
||||
{
|
||||
std::unordered_set<size_t> parts_ids;
|
||||
for (const auto &link : slice.overlaps_below) {
|
||||
size_t part_id = active_object_parts.get_flat_id(previous_slice_mappings.index_to_object_part_mapping.at(link.slice_idx));
|
||||
parts_ids.insert(part_id);
|
||||
transfered_weakest_connection.add(previous_slice_mappings.index_to_weakest_connection.at(link.slice_idx));
|
||||
}
|
||||
|
||||
final_part_id = *parts_ids.begin();
|
||||
for (size_t part_id : parts_ids) {
|
||||
if (final_part_id != part_id) {
|
||||
auto object_part = active_object_parts.access(part_id);
|
||||
if (auto object = to_partial_object(object_part)) {
|
||||
partial_objects.push_back(std::move(*object));
|
||||
}
|
||||
active_object_parts.merge(part_id, final_part_id);
|
||||
}
|
||||
}
|
||||
}
|
||||
const float bottom_z = layer->bottom_z();
|
||||
auto estimate_conn_strength = [bottom_z](const SliceConnection &conn) {
|
||||
if (conn.area < EPSILON) { // connection is empty, does not exists. Return max strength so that it is not picked as the
|
||||
// weakest connection.
|
||||
return INFINITY;
|
||||
}
|
||||
Vec3f centroid = conn.centroid_accumulator / conn.area;
|
||||
Vec2f variance = (conn.second_moment_of_area_accumulator / conn.area -
|
||||
centroid.head<2>().cwiseProduct(centroid.head<2>()));
|
||||
float xy_variance = variance.x() + variance.y();
|
||||
float arm_len_estimate = std::max(1.0f, bottom_z - (conn.centroid_accumulator.z() / conn.area));
|
||||
return conn.area * sqrt(xy_variance) / arm_len_estimate;
|
||||
};
|
||||
|
||||
#ifdef DETAILED_DEBUG_LOGS
|
||||
connection_to_below.print_info("new_weakest_connection");
|
||||
transfered_weakest_connection.print_info("transfered_weakest_connection");
|
||||
#endif
|
||||
|
||||
if (estimate_conn_strength(transfered_weakest_connection) > estimate_conn_strength(connection_to_below)) {
|
||||
transfered_weakest_connection = connection_to_below;
|
||||
}
|
||||
new_slice_mappings.index_to_weakest_connection.emplace(slice_idx, transfered_weakest_connection);
|
||||
new_slice_mappings.index_to_object_part_mapping.emplace(slice_idx, final_part_id);
|
||||
ObjectPart &part = active_object_parts.access(final_part_id);
|
||||
part.add(new_part);
|
||||
}
|
||||
}
|
||||
return new_slice_mappings;
|
||||
}
|
||||
|
||||
void reckon_global_supports(const tbb::concurrent_vector<ExtrusionLine> &external_perimeter_lines,
|
||||
const coordf_t layer_bottom_z,
|
||||
const Params ¶ms,
|
||||
ObjectPart &part,
|
||||
SliceConnection &weakest_connection,
|
||||
SupportPoints &supp_points,
|
||||
SupportGridFilter &supports_presence_grid)
|
||||
{
|
||||
LD current_slice_lines_distancer({external_perimeter_lines.begin(), external_perimeter_lines.end()});
|
||||
float unchecked_dist = params.min_distance_between_support_points + 1.0f;
|
||||
|
||||
for (const ExtrusionLine &line : external_perimeter_lines) {
|
||||
if ((unchecked_dist + line.len < params.min_distance_between_support_points &&
|
||||
line.curled_up_height < params.curling_tolerance_limit) ||
|
||||
line.len < EPSILON) {
|
||||
unchecked_dist += line.len;
|
||||
} else {
|
||||
unchecked_dist = line.len;
|
||||
Vec2f pivot_site_search_point = Vec2f(line.b + (line.b - line.a).normalized() * 300.0f);
|
||||
auto [dist, nidx, nearest_point] = current_slice_lines_distancer.distance_from_lines_extra<false>(pivot_site_search_point);
|
||||
Vec3f position = to_3d(nearest_point, layer_bottom_z);
|
||||
auto [force, cause] = part.is_stable_while_extruding(weakest_connection, line, position, layer_bottom_z, params);
|
||||
if (force > 0) {
|
||||
SupportPoint support_point{cause, position, params.support_points_interface_radius};
|
||||
reckon_new_support_point(part, weakest_connection, supp_points, supports_presence_grid, support_point, true);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::tuple<SupportPoints, PartialObjects> check_stability(const PrintObject *po,
|
||||
const PrecomputedSliceConnections &precomputed_slices_connections,
|
||||
const PrintTryCancel &cancel_func,
|
||||
@ -772,257 +1050,55 @@ std::tuple<SupportPoints, PartialObjects> check_stability(const PrintObject
|
||||
PartialObjects partial_objects{};
|
||||
LD prev_layer_ext_perim_lines;
|
||||
|
||||
std::unordered_map<size_t, size_t> prev_slice_idx_to_object_part_mapping;
|
||||
std::unordered_map<size_t, size_t> next_slice_idx_to_object_part_mapping;
|
||||
std::unordered_map<size_t, SliceConnection> prev_slice_idx_to_weakest_connection;
|
||||
std::unordered_map<size_t, SliceConnection> next_slice_idx_to_weakest_connection;
|
||||
SliceMappings slice_mappings;
|
||||
|
||||
auto remember_partial_object = [&active_object_parts, &partial_objects](size_t object_part_id) {
|
||||
auto object_part = active_object_parts.access(object_part_id);
|
||||
if (object_part.volume > EPSILON) {
|
||||
partial_objects.emplace_back(object_part.volume_centroid_accumulator / object_part.volume, object_part.volume,
|
||||
object_part.connected_to_bed);
|
||||
}
|
||||
};
|
||||
|
||||
for (size_t layer_idx = 0; layer_idx < po->layer_count(); ++layer_idx) {
|
||||
cancel_func();
|
||||
const Layer *layer = po->get_layer(layer_idx);
|
||||
float bottom_z = layer->bottom_z();
|
||||
auto create_support_point_position = [bottom_z](const Vec2f &layer_pos) { return Vec3f{layer_pos.x(), layer_pos.y(), bottom_z}; };
|
||||
|
||||
for (size_t slice_idx = 0; slice_idx < layer->lslices_ex.size(); ++slice_idx) {
|
||||
const LayerSlice &slice = layer->lslices_ex.at(slice_idx);
|
||||
const std::vector<const ExtrusionEntityCollection*> extrusion_collections{gather_extrusions(slice, layer)};
|
||||
const bool connected_to_bed = int(layer->id()) == params.raft_layers_count;
|
||||
slice_mappings = update_active_object_parts(layer, params, precomputed_slices_connections[layer_idx], slice_mappings, active_object_parts, partial_objects);
|
||||
|
||||
const std::optional<Polygons> brim{
|
||||
has_brim(layer, params) ?
|
||||
std::optional{get_brim(layer->lslices[slice_idx], params.brim_type, params.brim_width)} :
|
||||
std::nullopt
|
||||
};
|
||||
ObjectPart new_part{
|
||||
extrusion_collections,
|
||||
connected_to_bed,
|
||||
layer->print_z,
|
||||
layer->height,
|
||||
brim
|
||||
};
|
||||
std::optional<Linesf> prev_layer_boundary = layer->lower_layer != nullptr ?
|
||||
std::optional{to_unscaled_linesf(layer->lower_layer->lslices)} :
|
||||
std::nullopt;
|
||||
|
||||
const SliceConnection &connection_to_below = precomputed_slices_connections[layer_idx][slice_idx];
|
||||
|
||||
#ifdef DETAILED_DEBUG_LOGS
|
||||
std::cout << "SLICE IDX: " << slice_idx << std::endl;
|
||||
for (const auto &link : slice.overlaps_below) {
|
||||
std::cout << "connected to slice below: " << link.slice_idx << " by area : " << link.area << std::endl;
|
||||
}
|
||||
connection_to_below.print_info("CONNECTION TO BELOW");
|
||||
#endif
|
||||
|
||||
if (connection_to_below.area < EPSILON) { // new object part emerging
|
||||
size_t part_id = active_object_parts.insert(new_part);
|
||||
next_slice_idx_to_object_part_mapping.emplace(slice_idx, part_id);
|
||||
next_slice_idx_to_weakest_connection.emplace(slice_idx, connection_to_below);
|
||||
} else {
|
||||
size_t final_part_id{};
|
||||
SliceConnection transfered_weakest_connection{};
|
||||
// MERGE parts
|
||||
{
|
||||
std::unordered_set<size_t> parts_ids;
|
||||
for (const auto &link : slice.overlaps_below) {
|
||||
size_t part_id = active_object_parts.get_flat_id(prev_slice_idx_to_object_part_mapping.at(link.slice_idx));
|
||||
parts_ids.insert(part_id);
|
||||
transfered_weakest_connection.add(prev_slice_idx_to_weakest_connection.at(link.slice_idx));
|
||||
}
|
||||
|
||||
final_part_id = *parts_ids.begin();
|
||||
for (size_t part_id : parts_ids) {
|
||||
if (final_part_id != part_id) {
|
||||
remember_partial_object(part_id);
|
||||
active_object_parts.merge(part_id, final_part_id);
|
||||
}
|
||||
}
|
||||
}
|
||||
auto estimate_conn_strength = [bottom_z](const SliceConnection &conn) {
|
||||
if (conn.area < EPSILON) { // connection is empty, does not exists. Return max strength so that it is not picked as the
|
||||
// weakest connection.
|
||||
return INFINITY;
|
||||
}
|
||||
Vec3f centroid = conn.centroid_accumulator / conn.area;
|
||||
Vec2f variance = (conn.second_moment_of_area_accumulator / conn.area -
|
||||
centroid.head<2>().cwiseProduct(centroid.head<2>()));
|
||||
float xy_variance = variance.x() + variance.y();
|
||||
float arm_len_estimate = std::max(1.0f, bottom_z - (conn.centroid_accumulator.z() / conn.area));
|
||||
return conn.area * sqrt(xy_variance) / arm_len_estimate;
|
||||
};
|
||||
|
||||
#ifdef DETAILED_DEBUG_LOGS
|
||||
connection_to_below.print_info("new_weakest_connection");
|
||||
transfered_weakest_connection.print_info("transfered_weakest_connection");
|
||||
#endif
|
||||
|
||||
if (estimate_conn_strength(transfered_weakest_connection) > estimate_conn_strength(connection_to_below)) {
|
||||
transfered_weakest_connection = connection_to_below;
|
||||
}
|
||||
next_slice_idx_to_weakest_connection.emplace(slice_idx, transfered_weakest_connection);
|
||||
next_slice_idx_to_object_part_mapping.emplace(slice_idx, final_part_id);
|
||||
ObjectPart &part = active_object_parts.access(final_part_id);
|
||||
part.add(new_part);
|
||||
}
|
||||
}
|
||||
|
||||
prev_slice_idx_to_object_part_mapping = next_slice_idx_to_object_part_mapping;
|
||||
next_slice_idx_to_object_part_mapping.clear();
|
||||
prev_slice_idx_to_weakest_connection = next_slice_idx_to_weakest_connection;
|
||||
next_slice_idx_to_weakest_connection.clear();
|
||||
|
||||
auto get_flat_entities = [](const ExtrusionEntity *e) {
|
||||
std::vector<const ExtrusionEntity *> entities;
|
||||
std::vector<const ExtrusionEntity *> queue{e};
|
||||
while (!queue.empty()) {
|
||||
const ExtrusionEntity *next = queue.back();
|
||||
queue.pop_back();
|
||||
if (next->is_collection()) {
|
||||
for (const ExtrusionEntity *e : static_cast<const ExtrusionEntityCollection *>(next)->entities) {
|
||||
queue.push_back(e);
|
||||
}
|
||||
} else {
|
||||
entities.push_back(next);
|
||||
}
|
||||
}
|
||||
return entities;
|
||||
};
|
||||
|
||||
struct EnitityToCheck
|
||||
{
|
||||
const ExtrusionEntity *e;
|
||||
const LayerRegion *region;
|
||||
size_t slice_idx;
|
||||
};
|
||||
std::vector<EnitityToCheck> entities_to_check;
|
||||
for (size_t slice_idx = 0; slice_idx < layer->lslices_ex.size(); ++slice_idx) {
|
||||
const LayerSlice &slice = layer->lslices_ex.at(slice_idx);
|
||||
for (const auto &island : slice.islands) {
|
||||
for (const LayerExtrusionRange &fill_range : island.fills) {
|
||||
const LayerRegion *fill_region = layer->get_region(fill_range.region());
|
||||
for (size_t fill_idx : fill_range) {
|
||||
for (const ExtrusionEntity *e : get_flat_entities(fill_region->fills().entities[fill_idx])) {
|
||||
if (e->role() == ExtrusionRole::BridgeInfill) {
|
||||
entities_to_check.push_back({e, fill_region, slice_idx});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const LayerRegion *perimeter_region = layer->get_region(island.perimeters.region());
|
||||
for (size_t perimeter_idx : island.perimeters) {
|
||||
for (const ExtrusionEntity *e : get_flat_entities(perimeter_region->perimeters().entities[perimeter_idx])) {
|
||||
entities_to_check.push_back({e, perimeter_region, slice_idx});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
AABBTreeLines::LinesDistancer<Linef> prev_layer_boundary = layer->lower_layer != nullptr ?
|
||||
AABBTreeLines::LinesDistancer<Linef>{
|
||||
to_unscaled_linesf(layer->lower_layer->lslices)} :
|
||||
AABBTreeLines::LinesDistancer<Linef>{};
|
||||
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> unstable_lines_per_slice(layer->lslices_ex.size());
|
||||
std::vector<tbb::concurrent_vector<ExtrusionLine>> ext_perim_lines_per_slice(layer->lslices_ex.size());
|
||||
|
||||
tbb::parallel_for(tbb::blocked_range<size_t>(0, entities_to_check.size()),
|
||||
[&entities_to_check, &prev_layer_ext_perim_lines, &prev_layer_boundary, &unstable_lines_per_slice,
|
||||
&ext_perim_lines_per_slice, ¶ms](tbb::blocked_range<size_t> r) {
|
||||
for (size_t entity_idx = r.begin(); entity_idx < r.end(); ++entity_idx) {
|
||||
const auto &e_to_check = entities_to_check[entity_idx];
|
||||
for (const auto &line :
|
||||
check_extrusion_entity_stability(e_to_check.e, e_to_check.region, prev_layer_ext_perim_lines,
|
||||
prev_layer_boundary, params)) {
|
||||
if (line.support_point_generated.has_value()) {
|
||||
unstable_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
if (line.is_external_perimeter()) {
|
||||
ext_perim_lines_per_slice[e_to_check.slice_idx].push_back(line);
|
||||
}
|
||||
}
|
||||
}
|
||||
});
|
||||
LocalSupports local_supports{
|
||||
compute_local_supports(gather_entities_to_check(layer), prev_layer_boundary, prev_layer_ext_perim_lines, layer->lslices_ex.size(), params)};
|
||||
|
||||
std::vector<ExtrusionLine> current_layer_ext_perims_lines{};
|
||||
current_layer_ext_perims_lines.reserve(prev_layer_ext_perim_lines.get_lines().size());
|
||||
// All object parts updated, and for each slice we have coresponding weakest connection.
|
||||
// We can now check each slice and its corresponding weakest connection and object part for stability.
|
||||
for (size_t slice_idx = 0; slice_idx < layer->lslices_ex.size(); ++slice_idx) {
|
||||
ObjectPart &part = active_object_parts.access(prev_slice_idx_to_object_part_mapping[slice_idx]);
|
||||
SliceConnection &weakest_conn = prev_slice_idx_to_weakest_connection[slice_idx];
|
||||
ObjectPart &part = active_object_parts.access(slice_mappings.index_to_object_part_mapping[slice_idx]);
|
||||
SliceConnection &weakest_conn = slice_mappings.index_to_weakest_connection[slice_idx];
|
||||
|
||||
#ifdef DETAILED_DEBUG_LOGS
|
||||
weakest_conn.print_info("weakest connection info: ");
|
||||
#endif
|
||||
// Function that is used when new support point is generated. It will update the ObjectPart stability, weakest conneciton info,
|
||||
// and the support presence grid and add the point to the issues.
|
||||
auto reckon_new_support_point = [&part, &weakest_conn, &supp_points, &supports_presence_grid, ¶ms,
|
||||
&layer_idx](SupportPointCause cause, const Vec3f &support_point, float force,
|
||||
const Vec2f &dir) {
|
||||
// if position is taken and point is for global stability (force > 0) or we are too close to the bed, do not add
|
||||
// This allows local support points (e.g. bridging) to be generated densely
|
||||
if ((supports_presence_grid.position_taken(support_point) && force > 0) || layer_idx <= 1) {
|
||||
return;
|
||||
}
|
||||
|
||||
float area = params.support_points_interface_radius * params.support_points_interface_radius * float(PI);
|
||||
// add the stability effect of the point only if the spot is not taken, so that the densely created local support points do
|
||||
// not add unrealistic amount of stability to the object (due to overlaping of local support points)
|
||||
if (!(supports_presence_grid.position_taken(support_point))) {
|
||||
part.add_support_point(support_point, area);
|
||||
}
|
||||
|
||||
float radius = params.support_points_interface_radius;
|
||||
supp_points.emplace_back(cause, support_point, force, radius, dir);
|
||||
supports_presence_grid.take_position(support_point);
|
||||
|
||||
// The support point also increases the stability of the weakest connection of the object, which should be reflected
|
||||
if (weakest_conn.area > EPSILON) { // Do not add it to the weakest connection if it is not valid - does not exist
|
||||
weakest_conn.area += area;
|
||||
weakest_conn.centroid_accumulator += support_point * area;
|
||||
weakest_conn.second_moment_of_area_accumulator += area * support_point.head<2>().cwiseProduct(support_point.head<2>());
|
||||
weakest_conn.second_moment_of_area_covariance_accumulator += area * support_point.x() * support_point.y();
|
||||
}
|
||||
};
|
||||
|
||||
for (const auto &l : unstable_lines_per_slice[slice_idx]) {
|
||||
assert(l.support_point_generated.has_value());
|
||||
reckon_new_support_point(*l.support_point_generated, create_support_point_position(l.b), float(-EPSILON), Vec2f::Zero());
|
||||
}
|
||||
|
||||
LD current_slice_lines_distancer({ext_perim_lines_per_slice[slice_idx].begin(), ext_perim_lines_per_slice[slice_idx].end()});
|
||||
float unchecked_dist = params.min_distance_between_support_points + 1.0f;
|
||||
|
||||
for (const ExtrusionLine &line : current_slice_lines_distancer.get_lines()) {
|
||||
if ((unchecked_dist + line.len < params.min_distance_between_support_points && line.curled_up_height < params.curling_tolerance_limit) ||
|
||||
line.len < EPSILON) {
|
||||
unchecked_dist += line.len;
|
||||
} else {
|
||||
unchecked_dist = line.len;
|
||||
Vec2f pivot_site_search_point = Vec2f(line.b + (line.b - line.a).normalized() * 300.0f);
|
||||
auto [dist, nidx,
|
||||
nearest_point] = current_slice_lines_distancer.distance_from_lines_extra<false>(pivot_site_search_point);
|
||||
Vec3f support_point = create_support_point_position(nearest_point);
|
||||
auto [force, cause] = part.is_stable_while_extruding(weakest_conn, line, support_point, bottom_z, params);
|
||||
if (force > 0) {
|
||||
reckon_new_support_point(cause, support_point, force, (line.b - line.a).normalized());
|
||||
}
|
||||
if (layer_idx > 1) {
|
||||
for (const auto &l : local_supports.unstable_lines_per_slice[slice_idx]) {
|
||||
assert(l.support_point_generated.has_value());
|
||||
SupportPoint support_point{*l.support_point_generated, to_3d(l.b, bottom_z),
|
||||
params.support_points_interface_radius};
|
||||
reckon_new_support_point(part, weakest_conn, supp_points, supports_presence_grid, support_point);
|
||||
}
|
||||
}
|
||||
current_layer_ext_perims_lines.insert(current_layer_ext_perims_lines.end(), current_slice_lines_distancer.get_lines().begin(),
|
||||
current_slice_lines_distancer.get_lines().end());
|
||||
|
||||
const tbb::concurrent_vector<ExtrusionLine> &external_perimeter_lines = local_supports.ext_perim_lines_per_slice[slice_idx];
|
||||
if (layer_idx > 1) {
|
||||
reckon_global_supports(external_perimeter_lines, bottom_z, params, part, weakest_conn, supp_points, supports_presence_grid);
|
||||
}
|
||||
|
||||
current_layer_ext_perims_lines.insert(current_layer_ext_perims_lines.end(), external_perimeter_lines.begin(), external_perimeter_lines.end());
|
||||
} // slice iterations
|
||||
prev_layer_ext_perim_lines = LD(current_layer_ext_perims_lines);
|
||||
} // layer iterations
|
||||
|
||||
for (const auto& active_obj_pair : prev_slice_idx_to_object_part_mapping) {
|
||||
remember_partial_object(active_obj_pair.second);
|
||||
for (const auto& active_obj_pair : slice_mappings.index_to_object_part_mapping) {
|
||||
auto object_part = active_object_parts.access(active_obj_pair.second);
|
||||
if (auto object = to_partial_object(object_part)) {
|
||||
partial_objects.push_back(std::move(*object));
|
||||
}
|
||||
}
|
||||
|
||||
return {supp_points, partial_objects};
|
||||
|
||||
@ -104,8 +104,8 @@ enum class SupportPointCause {
|
||||
// Note that the force is only the difference - the amount needed to stabilize the object again.
|
||||
struct SupportPoint
|
||||
{
|
||||
SupportPoint(SupportPointCause cause, const Vec3f &position, float force, float spot_radius, const Vec2f &direction)
|
||||
: cause(cause), position(position), force(force), spot_radius(spot_radius), direction(direction)
|
||||
SupportPoint(SupportPointCause cause, const Vec3f &position, float spot_radius)
|
||||
: cause(cause), position(position), spot_radius(spot_radius)
|
||||
{}
|
||||
|
||||
bool is_local_extrusion_support() const
|
||||
@ -117,17 +117,9 @@ struct SupportPoint
|
||||
SupportPointCause cause; // reason why this support point was generated. Used for the user alerts
|
||||
// position is in unscaled coords. The z coordinate is aligned with the layers bottom_z coordiantes
|
||||
Vec3f position;
|
||||
// force that destabilizes the object to the point of falling/breaking. g*mm/s^2 units
|
||||
// It is valid only for global_object_support. For local extrusion support points, the force is -EPSILON
|
||||
// values gathered from large XL model: Min : 0 | Max : 18713800 | Average : 1361186 | Median : 329103
|
||||
// For reference 18713800 is weight of 1.8 Kg object, 329103 is weight of 0.03 Kg
|
||||
// The final sliced object weight was approx 0.5 Kg
|
||||
float force;
|
||||
// Expected spot size. The support point strength is calculated from the area defined by this value.
|
||||
// Currently equal to the support_points_interface_radius parameter above
|
||||
float spot_radius;
|
||||
// direction of the fall of the object (z part is neglected)
|
||||
Vec2f direction;
|
||||
};
|
||||
|
||||
using SupportPoints = std::vector<SupportPoint>;
|
||||
|
||||
@ -202,10 +202,10 @@ static void check_nfp(const std::string & outfile_prefix,
|
||||
ExPolygons bed_negative = diff_ex(bedrect, bedpoly);
|
||||
ExPolygons orb_ex_r = to_expolygons(orbiter);
|
||||
ExPolygons orb_ex_r_ch = {ExPolygon(Geometry::convex_hull(orb_ex_r))};
|
||||
auto orb_ex_offs_pos_r = offset_ex(orb_ex_r, SCALED_EPSILON);
|
||||
auto orb_ex_offs_neg_r = offset_ex(orb_ex_r, -SCALED_EPSILON);
|
||||
auto orb_ex_offs_pos_r_ch = offset_ex(orb_ex_r_ch, SCALED_EPSILON);
|
||||
auto orb_ex_offs_neg_r_ch = offset_ex(orb_ex_r_ch, -SCALED_EPSILON);
|
||||
auto orb_ex_offs_pos_r = offset_ex(orb_ex_r, scaled<float>(EPSILON));
|
||||
auto orb_ex_offs_neg_r = offset_ex(orb_ex_r, -scaled<float>(EPSILON));
|
||||
auto orb_ex_offs_pos_r_ch = offset_ex(orb_ex_r_ch, scaled<float>(EPSILON));
|
||||
auto orb_ex_offs_neg_r_ch = offset_ex(orb_ex_r_ch, -scaled<float>(EPSILON));
|
||||
|
||||
auto bedpoly_offs = offset_ex(bedpoly, SCALED_EPSILON);
|
||||
|
||||
|
||||
@ -738,8 +738,8 @@ bool is_collision_free(const Slic3r::Range<It> &item_range)
|
||||
|
||||
bool collision_free = true;
|
||||
foreach_combo(item_range, [&collision_free](auto &itm1, auto &itm2) {
|
||||
auto outline1 = offset(arr2::fixed_outline(itm1), -SCALED_EPSILON);
|
||||
auto outline2 = offset(arr2::fixed_outline(itm2), -SCALED_EPSILON);
|
||||
auto outline1 = offset(arr2::fixed_outline(itm1), -scaled<float>(EPSILON));
|
||||
auto outline2 = offset(arr2::fixed_outline(itm2), -scaled<float>(EPSILON));
|
||||
|
||||
auto inters = intersection(outline1, outline2);
|
||||
collision_free = collision_free && inters.empty();
|
||||
|
||||
@ -42,6 +42,7 @@ add_executable(${_TEST_NAME}_tests
|
||||
test_support_spots_generator.cpp
|
||||
../data/prusaparts.cpp
|
||||
../data/prusaparts.hpp
|
||||
test_static_map.cpp
|
||||
)
|
||||
|
||||
if (TARGET OpenVDB::openvdb)
|
||||
|
||||
@ -1,7 +1,11 @@
|
||||
#include <catch2/catch.hpp>
|
||||
#include <test_utils.hpp>
|
||||
|
||||
#include <random>
|
||||
|
||||
#include <libslic3r/Geometry/ArcWelder.hpp>
|
||||
#include <libslic3r/Geometry/Circle.hpp>
|
||||
#include <libslic3r/SVG.hpp>
|
||||
#include <libslic3r/libslic3r.h>
|
||||
|
||||
using namespace Slic3r;
|
||||
@ -13,25 +17,34 @@ TEST_CASE("arc basics", "[ArcWelder]") {
|
||||
Vec2f p1{ 2000.f, 1000.f };
|
||||
Vec2f p2{ 1000.f, 2000.f };
|
||||
float r{ 1000.f };
|
||||
const double s = 1000.f / sqrt(2.);
|
||||
THEN("90 degrees arc, CCW") {
|
||||
Vec2f c = ArcWelder::arc_center(p1, p2, r, true);
|
||||
Vec2f m = ArcWelder::arc_middle_point(p1, p2, r, true);
|
||||
REQUIRE(is_approx(c, Vec2f{ 1000.f, 1000.f }));
|
||||
REQUIRE(ArcWelder::arc_angle(p1, p2, r) == Approx(0.5 * M_PI));
|
||||
REQUIRE(ArcWelder::arc_length(p1, p2, r) == Approx(r * 0.5 * M_PI).epsilon(0.001));
|
||||
REQUIRE(is_approx(m, Vec2f{ 1000.f + s, 1000.f + s }, 0.001f));
|
||||
}
|
||||
THEN("90 degrees arc, CW") {
|
||||
Vec2f c = ArcWelder::arc_center(p1, p2, r, false);
|
||||
Vec2f m = ArcWelder::arc_middle_point(p1, p2, r, false);
|
||||
REQUIRE(is_approx(c, Vec2f{ 2000.f, 2000.f }));
|
||||
REQUIRE(is_approx(m, Vec2f{ 2000.f - s, 2000.f - s }, 0.001f));
|
||||
}
|
||||
THEN("270 degrees arc, CCW") {
|
||||
Vec2f c = ArcWelder::arc_center(p1, p2, - r, true);
|
||||
Vec2f m = ArcWelder::arc_middle_point(p1, p2, - r, true);
|
||||
REQUIRE(is_approx(c, Vec2f{ 2000.f, 2000.f }));
|
||||
REQUIRE(ArcWelder::arc_angle(p1, p2, - r) == Approx(1.5 * M_PI));
|
||||
REQUIRE(ArcWelder::arc_length(p1, p2, - r) == Approx(r * 1.5 * M_PI).epsilon(0.001));
|
||||
REQUIRE(is_approx(m, Vec2f{ 2000.f + s, 2000.f + s }, 0.001f));
|
||||
}
|
||||
THEN("270 degrees arc, CW") {
|
||||
Vec2f c = ArcWelder::arc_center(p1, p2, - r, false);
|
||||
Vec2f m = ArcWelder::arc_middle_point(p1, p2, - r, false);
|
||||
REQUIRE(is_approx(c, Vec2f{ 1000.f, 1000.f }));
|
||||
REQUIRE(is_approx(m, Vec2f{ 1000.f - s, 1000.f - s }, 0.001f));
|
||||
}
|
||||
}
|
||||
WHEN("arc from { 1707.11f, 1707.11f } to { 1000.f, 2000.f }") {
|
||||
@ -126,6 +139,28 @@ TEST_CASE("arc discretization", "[ArcWelder]") {
|
||||
}
|
||||
}
|
||||
|
||||
void test_arc_fit_variance(const Point &p1, const Point &p2, const float r, const float r_fit, const bool ccw, const Points::const_iterator begin, const Points::const_iterator end)
|
||||
{
|
||||
using namespace Slic3r::Geometry;
|
||||
double variance = ArcWelder::arc_fit_variance(p1, p2, r, ccw, begin, end);
|
||||
double variance_fit = ArcWelder::arc_fit_variance(p1, p2, r_fit, ccw, begin, end);
|
||||
REQUIRE(variance_fit < variance);
|
||||
};
|
||||
|
||||
void test_arc_fit_max_deviation(const Point &p1, const Point &p2, const float r, const float r_fit, const bool ccw, const Points::const_iterator begin, const Points::const_iterator end)
|
||||
{
|
||||
using namespace Slic3r::Geometry;
|
||||
double max_deviation = ArcWelder::arc_fit_max_deviation(p1, p2, r, ccw, begin, end);
|
||||
double max_deviation_fit = ArcWelder::arc_fit_max_deviation(p1, p2, r_fit, ccw, begin, end);
|
||||
REQUIRE(std::abs(max_deviation_fit) < std::abs(max_deviation));
|
||||
};
|
||||
|
||||
void test_arc_fit(const Point &p1, const Point &p2, const float r, const float r_fit, const bool ccw, const Points::const_iterator begin, const Points::const_iterator end)
|
||||
{
|
||||
test_arc_fit_variance(p1, p2, r, r_fit, ccw, begin, end);
|
||||
test_arc_fit_max_deviation(p1, p2, r, r_fit, ccw, begin, end);
|
||||
};
|
||||
|
||||
TEST_CASE("arc fitting", "[ArcWelder]") {
|
||||
using namespace Slic3r::Geometry;
|
||||
|
||||
@ -142,8 +177,8 @@ TEST_CASE("arc fitting", "[ArcWelder]") {
|
||||
REQUIRE(path.front().point == p1);
|
||||
REQUIRE(path.front().radius == 0.f);
|
||||
REQUIRE(path.back().point == p2);
|
||||
REQUIRE(path.back().radius == Approx(r));
|
||||
REQUIRE(path.back().ccw() == ccw);
|
||||
test_arc_fit(p1, p2, r, path.back().radius, true, pts.begin(), pts.end());
|
||||
};
|
||||
THEN("90 degrees arc, CCW is fitted") {
|
||||
test(p1, p2, radius, true);
|
||||
@ -169,6 +204,7 @@ TEST_CASE("arc fitting", "[ArcWelder]") {
|
||||
const float resolution = scaled<float>(0.002);
|
||||
auto test = [center1, center2, resolution](const Point &p1, const Point &p2, const Point &p3, const float r, const bool ccw) {
|
||||
Points pts = ArcWelder::arc_discretize(p1, p2, r, ccw, resolution);
|
||||
size_t num_pts1 = pts.size();
|
||||
{
|
||||
Points pts2 = ArcWelder::arc_discretize(p2, p3, - r, ! ccw, resolution);
|
||||
REQUIRE(pts.back() == pts2.front());
|
||||
@ -179,11 +215,11 @@ TEST_CASE("arc fitting", "[ArcWelder]") {
|
||||
REQUIRE(path.front().point == p1);
|
||||
REQUIRE(path.front().radius == 0.f);
|
||||
REQUIRE(path[1].point == p2);
|
||||
REQUIRE(path[1].radius == Approx(r));
|
||||
REQUIRE(path[1].ccw() == ccw);
|
||||
REQUIRE(path.back().point == p3);
|
||||
REQUIRE(path.back().radius == Approx(- r));
|
||||
REQUIRE(path.back().ccw() == ! ccw);
|
||||
test_arc_fit(p1, p2, r, path[1].radius, ccw, pts.begin(), pts.begin() + num_pts1);
|
||||
test_arc_fit(p2, p3, - r, path.back().radius, ! ccw, pts.begin() + num_pts1 - 1, pts.end());
|
||||
};
|
||||
THEN("90 degrees arches, CCW are fitted") {
|
||||
test(p1, p2, p3, radius, true);
|
||||
@ -200,6 +236,106 @@ TEST_CASE("arc fitting", "[ArcWelder]") {
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("least squares arc fitting, interpolating end points", "[ArcWelder]") {
|
||||
using namespace Slic3r::Geometry;
|
||||
|
||||
// Generate bunch of random arches.
|
||||
const coord_t max_coordinate = scaled<coord_t>(sqrt(250. - 1.));
|
||||
static constexpr const double min_radius = scaled<double>(0.01);
|
||||
static constexpr const double max_radius = scaled<double>(250.);
|
||||
// static constexpr const double deviation = scaled<double>(0.5);
|
||||
static constexpr const double deviation = scaled<double>(0.1);
|
||||
// Seeded with a fixed seed, to be repeatable.
|
||||
std::mt19937 rng(867092346);
|
||||
std::uniform_int_distribution<int32_t> coord_sampler(0, int32_t(max_coordinate));
|
||||
std::uniform_real_distribution<double> angle_sampler(0.001, 2. * M_PI - 0.001);
|
||||
std::uniform_real_distribution<double> radius_sampler(min_radius, max_radius);
|
||||
std::uniform_int_distribution<int> num_samples_sampler(1, 100);
|
||||
auto test_arc_fitting = [&rng, &coord_sampler, &num_samples_sampler, &angle_sampler, &radius_sampler]() {
|
||||
auto sample_point = [&rng, &coord_sampler]() {
|
||||
return Vec2d(coord_sampler(rng), coord_sampler(rng));
|
||||
};
|
||||
// Start and end point of the arc:
|
||||
Vec2d center_pos = sample_point();
|
||||
double angle0 = angle_sampler(rng);
|
||||
double angle = angle_sampler(rng);
|
||||
double radius = radius_sampler(rng);
|
||||
Vec2d v1 = Eigen::Rotation2D(angle0) * Vec2d(1., 0.);
|
||||
Vec2d v2 = Eigen::Rotation2D(angle0 + angle) * Vec2d(1., 0.);
|
||||
Vec2d start_pos = center_pos + radius * v1;
|
||||
Vec2d end_pos = center_pos + radius * v2;
|
||||
std::vector<Vec2d> samples;
|
||||
size_t num_samples = num_samples_sampler(rng);
|
||||
for (size_t i = 0; i < num_samples; ++ i) {
|
||||
double sample_r = sqrt(std::uniform_real_distribution<double>(sqr(std::max(0., radius - deviation)), sqr(radius + deviation))(rng));
|
||||
double sample_a = std::uniform_real_distribution<double>(0., angle)(rng);
|
||||
Vec2d pt = center_pos + Eigen::Rotation2D(angle0 + sample_a) * Vec2d(sample_r, 0.);
|
||||
samples.emplace_back(pt);
|
||||
assert((pt - center_pos).norm() > radius - deviation - SCALED_EPSILON);
|
||||
assert((pt - center_pos).norm() < radius + deviation + SCALED_EPSILON);
|
||||
}
|
||||
// Vec2d new_center = ArcWelder::arc_fit_center_algebraic_ls(start_pos, end_pos, center_pos, samples.begin(), samples.end());
|
||||
THEN("Center is fitted correctly") {
|
||||
std::optional<Vec2d> new_center_opt = ArcWelder::arc_fit_center_gauss_newton_ls(start_pos, end_pos, center_pos, samples.begin(), samples.end(), 10);
|
||||
REQUIRE(new_center_opt);
|
||||
if (new_center_opt) {
|
||||
Vec2d new_center = *new_center_opt;
|
||||
double new_radius = (new_center - start_pos).norm();
|
||||
double total_deviation = 0;
|
||||
double new_total_deviation = 0;
|
||||
for (const Vec2d &s : samples) {
|
||||
total_deviation += sqr((s - center_pos).norm() - radius);
|
||||
new_total_deviation += sqr((s - new_center).norm() - radius);
|
||||
}
|
||||
total_deviation /= double(num_samples);
|
||||
new_total_deviation /= double(num_samples);
|
||||
REQUIRE(new_total_deviation <= total_deviation);
|
||||
|
||||
#if 0
|
||||
double dist = (center_pos - new_center).norm();
|
||||
printf("Radius: %lf, Angle: %lf deg, Samples: %d, Dist: %lf\n", unscaled<double>(radius), 180. * angle / M_PI, int(num_samples), unscaled<double>(dist));
|
||||
// REQUIRE(is_approx(center_pos, new_center, deviation));
|
||||
if (dist > scaled<double>(1.)) {
|
||||
static int irun = 0;
|
||||
char path[2048];
|
||||
sprintf(path, "d:\\temp\\debug\\circle-fit-%d.svg", irun++);
|
||||
Eigen::AlignedBox<double, 2> bbox(start_pos, end_pos);
|
||||
for (const Vec2d& sample : samples)
|
||||
bbox.extend(sample);
|
||||
bbox.extend(center_pos);
|
||||
Slic3r::SVG svg(path, BoundingBox(bbox.min().cast<coord_t>(), bbox.max().cast<coord_t>()).inflated(bbox.sizes().maxCoeff() * 0.05));
|
||||
Points arc_src;
|
||||
for (size_t i = 0; i <= 1000; ++i)
|
||||
arc_src.emplace_back((center_pos + Eigen::Rotation2D(angle0 + double(i) * angle / 1000.) * Vec2d(radius, 0.)).cast<coord_t>());
|
||||
svg.draw(Polyline(arc_src));
|
||||
Points arc_new;
|
||||
double new_arc_angle = ArcWelder::arc_angle(start_pos, end_pos, (new_center - start_pos).norm());
|
||||
for (size_t i = 0; i <= 1000; ++i)
|
||||
arc_new.emplace_back((new_center + Eigen::Rotation2D(double(i) * new_arc_angle / 1000.) * (start_pos - new_center)).cast<coord_t>());
|
||||
svg.draw(Polyline(arc_new), "magenta");
|
||||
svg.draw(Point(start_pos.cast<coord_t>()), "blue");
|
||||
svg.draw(Point(end_pos.cast<coord_t>()), "blue");
|
||||
svg.draw(Point(center_pos.cast<coord_t>()), "blue");
|
||||
for (const Vec2d &sample : samples)
|
||||
svg.draw(Point(sample.cast<coord_t>()), "red");
|
||||
svg.draw(Point(new_center.cast<coord_t>()), "magenta");
|
||||
}
|
||||
if (!is_approx(center_pos, new_center, scaled<double>(5.))) {
|
||||
printf("Failed\n");
|
||||
}
|
||||
#endif
|
||||
} else {
|
||||
printf("Fitting failed\n");
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
WHEN("Generating a random arc and randomized arc samples") {
|
||||
for (size_t i = 0; i < 1000; ++ i)
|
||||
test_arc_fitting();
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("arc wedge test", "[ArcWelder]") {
|
||||
using namespace Slic3r::Geometry;
|
||||
|
||||
@ -262,3 +398,114 @@ TEST_CASE("arc wedge test", "[ArcWelder]") {
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#if 0
|
||||
// For quantization
|
||||
//#include <libslic3r/GCode/GCodeWriter.hpp>
|
||||
|
||||
TEST_CASE("arc quantization", "[ArcWelder]") {
|
||||
using namespace Slic3r::Geometry;
|
||||
|
||||
WHEN("generating a bunch of random arches") {
|
||||
static constexpr const size_t len = 100000;
|
||||
static constexpr const coord_t max_coordinate = scaled<coord_t>(250.);
|
||||
static constexpr const float max_radius = scaled<float>(250.);
|
||||
ArcWelder::Segments path;
|
||||
path.reserve(len + 1);
|
||||
// Seeded with a fixed seed, to be repeatable.
|
||||
std::mt19937 rng(987432690);
|
||||
// Generate bunch of random arches.
|
||||
std::uniform_int_distribution<int32_t> coord_sampler(0, int32_t(max_coordinate));
|
||||
std::uniform_real_distribution<float> radius_sampler(- max_radius, max_radius);
|
||||
std::uniform_int_distribution<int> orientation_sampler(0, 1);
|
||||
path.push_back({ Point{coord_sampler(rng), coord_sampler(rng)}, 0, ArcWelder::Orientation::CCW });
|
||||
for (size_t i = 0; i < len; ++ i) {
|
||||
ArcWelder::Segment seg { Point{coord_sampler(rng), coord_sampler(rng)}, radius_sampler(rng), orientation_sampler(rng) ? ArcWelder::Orientation::CCW : ArcWelder::Orientation::CW };
|
||||
while ((seg.point.cast<double>() - path.back().point.cast<double>()).norm() > 2. * std::abs(seg.radius))
|
||||
seg = { Point{coord_sampler(rng), coord_sampler(rng)}, radius_sampler(rng), orientation_sampler(rng) ? ArcWelder::Orientation::CCW : ArcWelder::Orientation::CW };
|
||||
path.push_back(seg);
|
||||
}
|
||||
// Run the test, quantize coordinates and radius, find the maximum error of quantization comparing the two arches.
|
||||
struct ArcEval {
|
||||
double error;
|
||||
double radius;
|
||||
double angle;
|
||||
};
|
||||
std::vector<ArcEval> center_errors_R;
|
||||
center_errors_R.reserve(len);
|
||||
std::vector<double> center_errors_R_exact;
|
||||
center_errors_R_exact.reserve(len);
|
||||
std::vector<double> center_errors_IJ;
|
||||
center_errors_IJ.reserve(len);
|
||||
for (size_t i = 0; i < len; ++ i) {
|
||||
// Source arc:
|
||||
const Vec2d start_point = unscaled<double>(path[i].point);
|
||||
const Vec2d end_point = unscaled<double>(path[i + 1].point);
|
||||
const double radius = unscaled<double>(path[i + 1].radius);
|
||||
const bool ccw = path[i + 1].ccw();
|
||||
const Vec2d center = ArcWelder::arc_center(start_point, end_point, radius, ccw);
|
||||
{
|
||||
const double d1 = (start_point - center).norm();
|
||||
const double d2 = (end_point - center).norm();
|
||||
const double dx = (end_point - start_point).norm();
|
||||
assert(std::abs(d1 - std::abs(radius)) < EPSILON);
|
||||
assert(std::abs(d2 - std::abs(radius)) < EPSILON);
|
||||
}
|
||||
// Quantized arc:
|
||||
const Vec2d start_point_quantized { GCodeFormatter::quantize_xyzf(start_point.x()), GCodeFormatter::quantize_xyzf(start_point.y()) };
|
||||
const Vec2d end_point_quantized { GCodeFormatter::quantize_xyzf(end_point .x()), GCodeFormatter::quantize_xyzf(end_point .y()) };
|
||||
const double radius_quantized { GCodeFormatter::quantize_xyzf(radius) };
|
||||
const Vec2d center_quantized { GCodeFormatter::quantize_xyzf(center .x()), GCodeFormatter::quantize_xyzf(center .y()) };
|
||||
// Evaulate maximum error for the quantized arc given by the end points and radius.
|
||||
const Vec2d center_from_quantized = ArcWelder::arc_center(start_point_quantized, end_point_quantized, radius, ccw);
|
||||
const Vec2d center_reconstructed = ArcWelder::arc_center(start_point_quantized, end_point_quantized, radius_quantized, ccw);
|
||||
#if 0
|
||||
center_errors_R.push_back({
|
||||
std::abs((center_reconstructed - center).norm()),
|
||||
radius,
|
||||
ArcWelder::arc_angle(start_point, end_point, radius) * 180. / M_PI
|
||||
});
|
||||
if (center_errors_R.back().error > 0.15)
|
||||
printf("Fuj\n");
|
||||
#endif
|
||||
center_errors_R_exact.emplace_back(std::abs((center_from_quantized - center).norm()));
|
||||
if (center_errors_R_exact.back() > 0.15)
|
||||
printf("Fuj\n");
|
||||
center_errors_IJ.emplace_back(std::abs((center_quantized - center).norm()));
|
||||
if (center_errors_IJ.back() > 0.15)
|
||||
printf("Fuj\n");
|
||||
|
||||
// Adjust center of the arc to the quantized end points.
|
||||
Vec2d third_point = ArcWelder::arc_middle_point(start_point, end_point, radius, ccw);
|
||||
double third_point_radius = (third_point - center).norm();
|
||||
assert(std::abs(third_point_radius - std::abs(radius)) < EPSILON);
|
||||
std::optional<Vec2d> center_adjusted = try_circle_center(start_point_quantized, end_point_quantized, third_point, EPSILON);
|
||||
//assert(center_adjusted);
|
||||
if (center_adjusted) {
|
||||
const double radius_adjusted = (center_adjusted.value() - start_point_quantized).norm() * (radius > 0 ? 1. : -1.);
|
||||
const double radius_adjusted_quantized = GCodeFormatter::quantize_xyzf(radius_adjusted);
|
||||
// Evaulate maximum error for the quantized arc given by the end points and radius.
|
||||
const Vec2f center_reconstructed = ArcWelder::arc_center(start_point_quantized.cast<float>(), end_point_quantized.cast<float>(), float(radius_adjusted_quantized), ccw);
|
||||
double rtest = std::abs(radius_adjusted_quantized);
|
||||
double d1 = std::abs((center_reconstructed - start_point.cast<float>()).norm() - rtest);
|
||||
double d2 = std::abs((center_reconstructed - end_point.cast<float>()).norm() - rtest);
|
||||
double d3 = std::abs((center_reconstructed - third_point.cast<float>()).norm() - rtest);
|
||||
double d = std::max(d1, std::max(d2, d3));
|
||||
center_errors_R.push_back({
|
||||
d,
|
||||
radius,
|
||||
ArcWelder::arc_angle(start_point, end_point, radius) * 180. / M_PI
|
||||
});
|
||||
} else {
|
||||
printf("Adjusted circle is collinear.\n");
|
||||
}
|
||||
}
|
||||
std::sort(center_errors_R.begin(), center_errors_R.end(), [](auto l, auto r) { return l.error > r.error; });
|
||||
std::sort(center_errors_R_exact.begin(), center_errors_R_exact.end(), [](auto l, auto r) { return l > r; });
|
||||
std::sort(center_errors_IJ.begin(), center_errors_IJ.end(), [](auto l, auto r) { return l > r; });
|
||||
printf("Maximum error R: %lf\n", center_errors_R.back().error);
|
||||
printf("Maximum error R exact: %lf\n", center_errors_R_exact.back());
|
||||
printf("Maximum error IJ: %lf\n", center_errors_IJ.back());
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -323,6 +323,54 @@ SCENARIO("Circle Fit, TaubinFit with Newton's method", "[Geometry]") {
|
||||
}
|
||||
}
|
||||
|
||||
SCENARIO("Circle Fit, least squares by decomposition or by solving normal equation", "[Geometry]") {
|
||||
auto test_circle_fit = [](const Geometry::Circled &circle, const Vec2d ¢er, const double radius) {
|
||||
THEN("A center point matches.") {
|
||||
REQUIRE(is_approx(circle.center, center));
|
||||
}
|
||||
THEN("Radius matches") {
|
||||
REQUIRE(is_approx(circle.radius, radius));
|
||||
}
|
||||
};
|
||||
|
||||
GIVEN("A vector of Vec2ds arranged in a half-circle with approximately the same distance R from some point") {
|
||||
const Vec2d expected_center(-6., 0.);
|
||||
const double expected_radius = 6.;
|
||||
Vec2ds sample{Vec2d(6.0, 0), Vec2d(5.1961524, 3), Vec2d(3 ,5.1961524), Vec2d(0, 6.0), Vec2d(3, 5.1961524), Vec2d(-5.1961524, 3), Vec2d(-6.0, 0)};
|
||||
std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Vec2d &a) { return a + expected_center; });
|
||||
|
||||
WHEN("Circle fit is called on the entire array, least squares SVD") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_svd(sample), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the first four points, least squares SVD") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_svd(Vec2ds(sample.cbegin(), sample.cbegin() + 4)), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the middle four points, least squares SVD") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_svd(Vec2ds(sample.cbegin() + 2, sample.cbegin() + 6)), expected_center, expected_radius);
|
||||
}
|
||||
|
||||
WHEN("Circle fit is called on the entire array, least squares QR decomposition") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_qr(sample), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the first four points, least squares QR decomposition") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_qr(Vec2ds(sample.cbegin(), sample.cbegin() + 4)), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the middle four points, least squares QR decomposition") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_qr(Vec2ds(sample.cbegin() + 2, sample.cbegin() + 6)), expected_center, expected_radius);
|
||||
}
|
||||
|
||||
WHEN("Circle fit is called on the entire array, least squares by normal equations") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_normal(sample), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the first four points, least squares by normal equations") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_normal(Vec2ds(sample.cbegin(), sample.cbegin() + 4)), expected_center, expected_radius);
|
||||
}
|
||||
WHEN("Circle fit is called on the middle four points, least squares by normal equations") {
|
||||
test_circle_fit(Geometry::circle_linear_least_squares_normal(Vec2ds(sample.cbegin() + 2, sample.cbegin() + 6)), expected_center, expected_radius);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("smallest_enclosing_circle_welzl", "[Geometry]") {
|
||||
// Some random points in plane.
|
||||
Points pts {
|
||||
|
||||
98
tests/libslic3r/test_static_map.cpp
Normal file
98
tests/libslic3r/test_static_map.cpp
Normal file
@ -0,0 +1,98 @@
|
||||
#include <catch2/catch.hpp>
|
||||
#include <string_view>
|
||||
|
||||
#include "libslic3r/StaticMap.hpp"
|
||||
|
||||
TEST_CASE("Empty static map should be possible to create and should be empty", "[StaticMap]")
|
||||
{
|
||||
using namespace Slic3r;
|
||||
|
||||
static const constexpr StaticSet EmptySet;
|
||||
|
||||
static const constexpr auto EmptyMap = make_staticmap<int, int>();
|
||||
|
||||
constexpr bool is_map_empty = EmptyMap.empty();
|
||||
constexpr bool is_set_empty = EmptySet.empty();
|
||||
|
||||
REQUIRE(is_map_empty);
|
||||
REQUIRE(is_set_empty);
|
||||
}
|
||||
|
||||
TEST_CASE("StaticSet should derive it's type from the initializer", "[StaticMap]") {
|
||||
using namespace Slic3r;
|
||||
static const constexpr StaticSet iOneSet = { 1 };
|
||||
static constexpr size_t iOneSetSize = iOneSet.size();
|
||||
|
||||
REQUIRE(iOneSetSize == 1);
|
||||
|
||||
static const constexpr StaticSet iManySet = { 1, 3, 5, 80, 40 };
|
||||
static constexpr size_t iManySetSize = iManySet.size();
|
||||
|
||||
REQUIRE(iManySetSize == 5);
|
||||
}
|
||||
|
||||
TEST_CASE("StaticMap should derive it's type using make_staticmap", "[StaticMap]") {
|
||||
using namespace Slic3r;
|
||||
static const constexpr auto ciOneMap = make_staticmap<char, int>({
|
||||
{'a', 1},
|
||||
});
|
||||
|
||||
static constexpr size_t ciOneMapSize = ciOneMap.size();
|
||||
static constexpr bool ciOneMapValid = query(ciOneMap, 'a').value_or(0) == 1;
|
||||
|
||||
REQUIRE(ciOneMapSize == 1);
|
||||
REQUIRE(ciOneMapValid);
|
||||
|
||||
static const constexpr auto ciManyMap = make_staticmap<char, int>({
|
||||
{'a', 1}, {'b', 2}, {'A', 10}
|
||||
});
|
||||
|
||||
static constexpr size_t ciManyMapSize = ciManyMap.size();
|
||||
static constexpr bool ciManyMapValid =
|
||||
query(ciManyMap, 'a').value_or(0) == 1 &&
|
||||
query(ciManyMap, 'b').value_or(0) == 2 &&
|
||||
query(ciManyMap, 'A').value_or(0) == 10 &&
|
||||
!contains(ciManyMap, 'B') &&
|
||||
!query(ciManyMap, 'c').has_value();
|
||||
|
||||
REQUIRE(ciManyMapSize == 3);
|
||||
REQUIRE(ciManyMapValid);
|
||||
|
||||
for (auto &[k, v] : ciManyMap) {
|
||||
auto val = query(ciManyMap, k);
|
||||
REQUIRE(val.has_value());
|
||||
REQUIRE(*val == v);
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("StaticSet should be able to find contained values", "[StaticMap]")
|
||||
{
|
||||
using namespace Slic3r;
|
||||
using namespace std::string_view_literals;
|
||||
|
||||
auto cmp = [](const char *a, const char *b) constexpr {
|
||||
return std::string_view{a} < std::string_view{b};
|
||||
};
|
||||
|
||||
static constexpr StaticSet CStrSet = {cmp, "One", "Two", "Three"};
|
||||
static constexpr StaticSet StringSet = {"One"sv, "Two"sv, "Three"sv};
|
||||
|
||||
static constexpr bool CStrSetValid = query(CStrSet, "One").has_value() &&
|
||||
contains(CStrSet, "Two") &&
|
||||
contains(CStrSet, "Three") &&
|
||||
!contains(CStrSet, "one") &&
|
||||
!contains(CStrSet, "two") &&
|
||||
!contains(CStrSet, "three");
|
||||
|
||||
static constexpr bool StringSetValid = contains(StringSet, "One"sv) &&
|
||||
contains(StringSet, "Two"sv) &&
|
||||
contains(StringSet, "Three"sv) &&
|
||||
!contains(StringSet, "one"sv) &&
|
||||
!contains(StringSet, "two"sv) &&
|
||||
!contains(StringSet, "three"sv);
|
||||
|
||||
REQUIRE(CStrSetValid);
|
||||
REQUIRE(StringSetValid);
|
||||
REQUIRE(CStrSet.size() == 3);
|
||||
REQUIRE(StringSet.size() == 3);
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user