diff --git a/doc/seam_placement/corner_penalty_function.png b/doc/seam_placement/corner_penalty_function.png new file mode 100644 index 000000000..bf445d671 Binary files /dev/null and b/doc/seam_placement/corner_penalty_function.png differ diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp index 217166f8c..3ea18de8d 100644 --- a/src/libslic3r/AABBTreeIndirect.hpp +++ b/src/libslic3r/AABBTreeIndirect.hpp @@ -709,10 +709,15 @@ inline bool intersect_ray_all_hits( origin, dir, VectorType(dir.cwiseInverse()), eps } }; - if (! tree.empty()) { + if (tree.empty()) { + hits.clear(); + } else { + // Reusing the output memory if there is some memory already pre-allocated. + ray_intersector.hits = std::move(hits); + ray_intersector.hits.clear(); ray_intersector.hits.reserve(8); detail::intersect_ray_recursive_all_hits(ray_intersector, 0); - std::swap(hits, ray_intersector.hits); + hits = std::move(ray_intersector.hits); std::sort(hits.begin(), hits.end(), [](const auto &l, const auto &r) { return l.t < r.t; }); } return ! hits.empty(); @@ -759,8 +764,8 @@ inline bool is_any_triangle_in_radius( const TreeType &tree, // Point to which the closest point on the indexed triangle set is searched for. const VectorType &point, - // Maximum distance in which triangle is search for - typename VectorType::Scalar &max_distance) + //Square of maximum distance in which triangle is searched for + typename VectorType::Scalar &max_distance_squared) { using Scalar = typename VectorType::Scalar; auto distancer = detail::IndexedTriangleSetDistancer @@ -774,7 +779,7 @@ inline bool is_any_triangle_in_radius( return false; } - detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance, hit_idx, hit_point); + detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point); return hit_point.allFinite(); } diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 838af44fb..a9ce8dd1f 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -138,10 +138,12 @@ set(SLIC3R_SOURCES GCodeWriter.hpp Geometry.cpp Geometry.hpp + Geometry/Bicubic.hpp Geometry/Circle.cpp Geometry/Circle.hpp Geometry/ConvexHull.cpp Geometry/ConvexHull.hpp + Geometry/Curves.hpp Geometry/MedialAxis.cpp Geometry/MedialAxis.hpp Geometry/Voronoi.hpp @@ -225,6 +227,8 @@ set(SLIC3R_SOURCES SlicesToTriangleMesh.cpp SlicingAdaptive.cpp SlicingAdaptive.hpp + Subdivide.cpp + Subdivide.hpp SupportMaterial.cpp SupportMaterial.hpp Surface.cpp diff --git a/src/libslic3r/Color.hpp b/src/libslic3r/Color.hpp index 183705c4a..8044a0318 100644 --- a/src/libslic3r/Color.hpp +++ b/src/libslic3r/Color.hpp @@ -4,6 +4,8 @@ #include #include +#include "Point.hpp" + namespace Slic3r { class ColorRGB @@ -133,41 +135,57 @@ public: static const ColorRGBA Z() { return { 0.0f, 0.0f, 0.75f, 1.0f }; } }; -extern ColorRGB operator * (float value, const ColorRGB& other); -extern ColorRGBA operator * (float value, const ColorRGBA& other); +ColorRGB operator * (float value, const ColorRGB& other); +ColorRGBA operator * (float value, const ColorRGBA& other); -extern ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t); -extern ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t); +ColorRGB lerp(const ColorRGB& a, const ColorRGB& b, float t); +ColorRGBA lerp(const ColorRGBA& a, const ColorRGBA& b, float t); -extern ColorRGB complementary(const ColorRGB& color); -extern ColorRGBA complementary(const ColorRGBA& color); +ColorRGB complementary(const ColorRGB& color); +ColorRGBA complementary(const ColorRGBA& color); -extern ColorRGB saturate(const ColorRGB& color, float factor); -extern ColorRGBA saturate(const ColorRGBA& color, float factor); +ColorRGB saturate(const ColorRGB& color, float factor); +ColorRGBA saturate(const ColorRGBA& color, float factor); -extern ColorRGB opposite(const ColorRGB& color); -extern ColorRGB opposite(const ColorRGB& a, const ColorRGB& b); +ColorRGB opposite(const ColorRGB& color); +ColorRGB opposite(const ColorRGB& a, const ColorRGB& b); -extern bool can_decode_color(const std::string& color); +bool can_decode_color(const std::string& color); -extern bool decode_color(const std::string& color_in, ColorRGB& color_out); -extern bool decode_color(const std::string& color_in, ColorRGBA& color_out); +bool decode_color(const std::string& color_in, ColorRGB& color_out); +bool decode_color(const std::string& color_in, ColorRGBA& color_out); -extern bool decode_colors(const std::vector& colors_in, std::vector& colors_out); -extern bool decode_colors(const std::vector& colors_in, std::vector& colors_out); +bool decode_colors(const std::vector& colors_in, std::vector& colors_out); +bool decode_colors(const std::vector& colors_in, std::vector& colors_out); -extern std::string encode_color(const ColorRGB& color); -extern std::string encode_color(const ColorRGBA& color); +std::string encode_color(const ColorRGB& color); +std::string encode_color(const ColorRGBA& color); -extern ColorRGB to_rgb(const ColorRGBA& other_rgba); -extern ColorRGBA to_rgba(const ColorRGB& other_rgb); -extern ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha); +ColorRGB to_rgb(const ColorRGBA& other_rgba); +ColorRGBA to_rgba(const ColorRGB& other_rgb); +ColorRGBA to_rgba(const ColorRGB& other_rgb, float alpha); -extern ColorRGBA picking_decode(unsigned int id); -extern unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b); +// Color mapping of a value into RGB false colors. +inline Vec3f value_to_rgbf(float minimum, float maximum, float value) +{ + float ratio = 2.0f * (value - minimum) / (maximum - minimum); + float b = std::max(0.0f, (1.0f - ratio)); + float r = std::max(0.0f, (ratio - 1.0f)); + float g = 1.0f - b - r; + return Vec3f { r, g, b }; +} + +// Color mapping of a value into RGB false colors. +inline Vec3i value_to_rgbi(float minimum, float maximum, float value) +{ + return (value_to_rgbf(minimum, maximum, value) * 255).cast(); +} + +ColorRGBA picking_decode(unsigned int id); +unsigned int picking_encode(unsigned char r, unsigned char g, unsigned char b); // Produce an alpha channel checksum for the red green blue components. The alpha channel may then be used to verify, whether the rgb components // were not interpolated by alpha blending or multi sampling. -extern unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue); +unsigned char picking_checksum_alpha_channel(unsigned char red, unsigned char green, unsigned char blue); } // namespace Slic3r diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index 3d0ca5069..1a6ee4b80 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -2312,7 +2312,6 @@ GCode::LayerResult GCode::process_layer( } // for objects // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders. - std::vector> lower_layer_edge_grids(layers.size()); for (unsigned int extruder_id : layer_tools.extruders) { gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ? @@ -2337,7 +2336,7 @@ GCode::LayerResult GCode::process_layer( path.mm3_per_mm = mm3_per_mm; } //FIXME using the support_material_speed of the 1st object printed. - gcode += this->extrude_loop(loop, "skirt", m_config.support_material_speed.value); + gcode += this->extrude_loop(loop, "skirt"sv, m_config.support_material_speed.value); } m_avoid_crossing_perimeters.use_external_mp(false); // Allow a straight travel move to the first object point if this is the first layer (but don't in next layers). @@ -2350,7 +2349,7 @@ GCode::LayerResult GCode::process_layer( this->set_origin(0., 0.); m_avoid_crossing_perimeters.use_external_mp(); for (const ExtrusionEntity *ee : print.brim().entities) { - gcode += this->extrude_entity(*ee, "brim", m_config.support_material_speed.value); + gcode += this->extrude_entity(*ee, "brim"sv, m_config.support_material_speed.value); } m_brim_done = true; m_avoid_crossing_perimeters.use_external_mp(false); @@ -2406,9 +2405,9 @@ GCode::LayerResult GCode::process_layer( //FIXME the following code prints regions in the order they are defined, the path is not optimized in any way. if (print.config().infill_first) { gcode += this->extrude_infill(print, by_region_specific, false); - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); } else { - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); gcode += this->extrude_infill(print,by_region_specific, false); } // ironing @@ -2543,51 +2542,28 @@ std::string GCode::change_layer(coordf_t print_z) return gcode; } - - -static std::unique_ptr calculate_layer_edge_grid(const Layer& layer) -{ - auto out = make_unique(); - - // Create the distance field for a layer below. - const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5); - out->create(layer.lslices, distance_field_resolution); - out->calculate_sdf(); -#if 0 - { - static int iRun = 0; - BoundingBox bbox = (*lower_layer_edge_grid)->bbox(); - bbox.min(0) -= scale_(5.f); - bbox.min(1) -= scale_(5.f); - bbox.max(0) += scale_(5.f); - bbox.max(1) += scale_(5.f); - EdgeGrid::save_png(*(*lower_layer_edge_grid), bbox, scale_(0.1f), debug_out_path("GCode_extrude_loop_edge_grid-%d.png", iRun++)); - } -#endif - return out; +static const auto comment_perimeter = "perimeter"sv; +// Comparing string_view pointer & length for speed. +static inline bool comment_is_perimeter(const std::string_view comment) { + return comment.data() == comment_perimeter.data() && comment.size() == comment_perimeter.size(); } - -std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, double speed, std::unique_ptr *lower_layer_edge_grid) +std::string GCode::extrude_loop(ExtrusionLoop loop, const std::string_view description, double speed) { // get a copy; don't modify the orientation of the original loop object otherwise // next copies (if any) would not detect the correct orientation - if (m_layer->lower_layer && lower_layer_edge_grid != nullptr && ! *lower_layer_edge_grid) - *lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer); - // extrude all loops ccw bool was_clockwise = loop.make_counter_clockwise(); // find the point of the loop that is closest to the current extruder position // or randomize if requested Point last_pos = this->last_pos(); - if (m_config.spiral_vase) { + if (! m_config.spiral_vase && comment_is_perimeter(description)) { + assert(m_layer != nullptr); + m_seam_placer.place_seam(m_layer, loop, m_config.external_perimeters_first, this->last_pos()); + } else loop.split_at(last_pos, false); - } - else - m_seam_placer.place_seam(loop, this->last_pos(), m_config.external_perimeters_first, - EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr); // clip the path to avoid the extruder to get exactly on the first point of the loop; // if polyline was shorter than the clipping distance we'd get a null polyline, so @@ -2607,11 +2583,9 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou // extrude along the path std::string gcode; - for (ExtrusionPaths::iterator path = paths.begin(); path != paths.end(); ++path) { -// description += ExtrusionLoop::role_to_string(loop.loop_role()); -// description += ExtrusionEntity::role_to_string(path->role); - path->simplify(m_scaled_resolution); - gcode += this->_extrude(*path, description, speed); + for (ExtrusionPath &path : paths) { + path.simplify(m_scaled_resolution); + gcode += this->_extrude(path, description, speed); } // reset acceleration @@ -2626,18 +2600,19 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou // the side depends on the original winding order of the polygon (left for contours, right for holes) //FIXME improve the algorithm in case the loop is tiny. //FIXME improve the algorithm in case the loop is split into segments with a low number of points (see the Point b query). - Point a = paths.front().polyline.points[1]; // second point - Point b = *(paths.back().polyline.points.end()-3); // second to last point + // Angle from the 2nd point to the last point. + double angle_inside = angle(paths.front().polyline.points[1] - paths.front().first_point(), + *(paths.back().polyline.points.end()-3) - paths.front().first_point()); + assert(angle_inside >= -M_PI && angle_inside <= M_PI); + // 3rd of this angle will be taken, thus make the angle monotonic before interpolation. if (was_clockwise) { - // swap points - Point c = a; a = b; b = c; + if (angle_inside > 0) + angle_inside -= 2.0 * M_PI; + } else { + if (angle_inside < 0) + angle_inside += 2.0 * M_PI; } - double angle = paths.front().first_point().ccw_angle(a, b) / 3; - - // turn left if contour, turn right if hole - if (was_clockwise) angle *= -1; - // create the destination point along the first segment and rotate it // we make sure we don't exceed the segment length because we don't know // the rotation of the second segment so we might cross the object boundary @@ -2649,7 +2624,8 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou // Shift by no more than a nozzle diameter. //FIXME Hiding the seams will not work nicely for very densely discretized contours! Point pt = ((nd * nd >= l2) ? p2 : (p1 + v * (nd / sqrt(l2)))).cast(); - pt.rotate(angle, paths.front().polyline.points.front()); + // Rotate pt inside around the seam point. + pt.rotate(angle_inside / 3., paths.front().polyline.points.front()); // generate the travel move gcode += m_writer.travel_to_xy(this->point_to_gcode(pt), "move inwards before travel"); } @@ -2657,13 +2633,11 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou return gcode; } -std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string description, double speed) +std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed) { // extrude along the path std::string gcode; for (ExtrusionPath path : multipath.paths) { -// description += ExtrusionLoop::role_to_string(loop.loop_role()); -// description += ExtrusionEntity::role_to_string(path->role); path.simplify(m_scaled_resolution); gcode += this->_extrude(path, description, speed); } @@ -2676,22 +2650,21 @@ std::string GCode::extrude_multi_path(ExtrusionMultiPath multipath, std::string return gcode; } -std::string GCode::extrude_entity(const ExtrusionEntity &entity, std::string description, double speed, std::unique_ptr *lower_layer_edge_grid) +std::string GCode::extrude_entity(const ExtrusionEntity &entity, const std::string_view description, double speed) { if (const ExtrusionPath* path = dynamic_cast(&entity)) return this->extrude_path(*path, description, speed); else if (const ExtrusionMultiPath* multipath = dynamic_cast(&entity)) return this->extrude_multi_path(*multipath, description, speed); else if (const ExtrusionLoop* loop = dynamic_cast(&entity)) - return this->extrude_loop(*loop, description, speed, lower_layer_edge_grid); + return this->extrude_loop(*loop, description, speed); else throw Slic3r::InvalidArgument("Invalid argument supplied to extrude()"); return ""; } -std::string GCode::extrude_path(ExtrusionPath path, std::string description, double speed) +std::string GCode::extrude_path(ExtrusionPath path, std::string_view description, double speed) { -// description += ExtrusionEntity::role_to_string(path.role()); path.simplify(m_scaled_resolution); std::string gcode = this->_extrude(path, description, speed); if (m_wipe.enable) { @@ -2704,24 +2677,15 @@ std::string GCode::extrude_path(ExtrusionPath path, std::string description, dou } // Extrude perimeters: Decide where to put seams (hide or align seams). -std::string GCode::extrude_perimeters(const Print &print, const std::vector &by_region, std::unique_ptr &lower_layer_edge_grid) +std::string GCode::extrude_perimeters(const Print &print, const std::vector &by_region) { std::string gcode; for (const ObjectByExtruder::Island::Region ®ion : by_region) if (! region.perimeters.empty()) { m_config.apply(print.get_print_region(®ion - &by_region.front()).config()); - // plan_perimeters tries to place seams, it needs to have the lower_layer_edge_grid calculated already. - if (m_layer->lower_layer && ! lower_layer_edge_grid) - lower_layer_edge_grid = calculate_layer_edge_grid(*m_layer->lower_layer); - - m_seam_placer.plan_perimeters(std::vector(region.perimeters.begin(), region.perimeters.end()), - *m_layer, m_config.seam_position, this->last_pos(), EXTRUDER_CONFIG(nozzle_diameter), - (m_layer == NULL ? nullptr : m_layer->object()), - (lower_layer_edge_grid ? lower_layer_edge_grid.get() : nullptr)); - for (const ExtrusionEntity* ee : region.perimeters) - gcode += this->extrude_entity(*ee, "perimeter", -1., &lower_layer_edge_grid); + gcode += this->extrude_entity(*ee, comment_perimeter, -1.); } return gcode; } @@ -2731,7 +2695,7 @@ std::string GCode::extrude_infill(const Print &print, const std::vectorrole(); assert(role == erSupportMaterial || role == erSupportMaterialInterface); - const char *label = (role == erSupportMaterial) ? support_label : support_interface_label; + const auto label = (role == erSupportMaterial) ? support_label : support_interface_label; const double speed = (role == erSupportMaterial) ? support_speed : support_interface_speed; const ExtrusionPath *path = dynamic_cast(ee); if (path) @@ -2855,20 +2819,18 @@ void GCode::GCodeOutputStream::write_format(const char* format, ...) va_end(args); } -std::string GCode::_extrude(const ExtrusionPath &path, std::string description, double speed) +std::string GCode::_extrude(const ExtrusionPath &path, const std::string_view description, double speed) { std::string gcode; - - if (is_bridge(path.role())) - description += " (bridge)"; + const std::string_view description_bridge = is_bridge(path.role()) ? " (bridge)"sv : ""sv; // go to first point of extrusion path if (!m_last_pos_defined || m_last_pos != path.first_point()) { - gcode += this->travel_to( - path.first_point(), - path.role(), - "move to first " + description + " point" - ); + std::string comment = "move to first "; + comment += description; + comment += description_bridge; + comment += " point"; + gcode += this->travel_to(path.first_point(), path.role(), comment); } // compensate retraction @@ -3006,7 +2968,11 @@ std::string GCode::_extrude(const ExtrusionPath &path, std::string description, gcode += m_writer.set_speed(F, "", comment); double path_length = 0.; { - std::string comment = m_config.gcode_comments ? description : ""; + std::string comment; + if (m_config.gcode_comments) { + comment = description; + comment += description_bridge; + } Vec2d prev = this->point_to_gcode_quantized(path.polyline.points.front()); auto it = path.polyline.points.begin(); auto end = path.polyline.points.end(); diff --git a/src/libslic3r/GCode.hpp b/src/libslic3r/GCode.hpp index aa4682e34..5afe3d915 100644 --- a/src/libslic3r/GCode.hpp +++ b/src/libslic3r/GCode.hpp @@ -274,10 +274,10 @@ private: void set_extruders(const std::vector &extruder_ids); std::string preamble(); std::string change_layer(coordf_t print_z); - std::string extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1., std::unique_ptr *lower_layer_edge_grid = nullptr); - std::string extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1., std::unique_ptr *lower_layer_edge_grid = nullptr); - std::string extrude_multi_path(ExtrusionMultiPath multipath, std::string description = "", double speed = -1.); - std::string extrude_path(ExtrusionPath path, std::string description = "", double speed = -1.); + std::string extrude_entity(const ExtrusionEntity &entity, const std::string_view description, double speed = -1.); + std::string extrude_loop(ExtrusionLoop loop, const std::string_view description, double speed = -1.); + std::string extrude_multi_path(ExtrusionMultiPath multipath, const std::string_view description, double speed = -1.); + std::string extrude_path(ExtrusionPath path, const std::string_view description, double speed = -1.); // Extruding multiple objects with soluble / non-soluble / combined supports // on a multi-material printer, trying to minimize tool switches. @@ -342,7 +342,7 @@ private: // For sequential print, the instance of the object to be printing has to be defined. const size_t single_object_instance_idx); - std::string extrude_perimeters(const Print &print, const std::vector &by_region, std::unique_ptr &lower_layer_edge_grid); + std::string extrude_perimeters(const Print &print, const std::vector &by_region); std::string extrude_infill(const Print &print, const std::vector &by_region, bool ironing); std::string extrude_support(const ExtrusionEntityCollection &support_fills); @@ -428,7 +428,7 @@ private: // Processor GCodeProcessor m_processor; - std::string _extrude(const ExtrusionPath &path, std::string description = "", double speed = -1); + std::string _extrude(const ExtrusionPath &path, const std::string_view description, double speed = -1); void print_machine_envelope(GCodeOutputStream &file, Print &print); void _print_first_layer_bed_temperature(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait); void _print_first_layer_extruder_temperatures(GCodeOutputStream &file, Print &print, const std::string &gcode, unsigned int first_printing_extruder_id, bool wait); diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index bcd238a72..6e0968145 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -1,1011 +1,1416 @@ #include "SeamPlacer.hpp" +#include "tbb/parallel_for.h" +#include "tbb/blocked_range.h" +#include "tbb/parallel_reduce.h" +#include +#include +#include +#include + #include "libslic3r/ExtrusionEntity.hpp" #include "libslic3r/Print.hpp" #include "libslic3r/BoundingBox.hpp" +#include "libslic3r/Color.hpp" #include "libslic3r/EdgeGrid.hpp" #include "libslic3r/ClipperUtils.hpp" -#include "libslic3r/SVG.hpp" #include "libslic3r/Layer.hpp" +#include "libslic3r/QuadricEdgeCollapse.hpp" +#include "libslic3r/Subdivide.hpp" + +#include "libslic3r/Geometry/Curves.hpp" + +#include "libslic3r/Utils.hpp" + +//#define DEBUG_FILES + +#ifdef DEBUG_FILES +#include +#include +#endif namespace Slic3r { -// This penalty is added to all points inside custom blockers (subtracted from pts inside enforcers). -static constexpr float ENFORCER_BLOCKER_PENALTY = 100; +namespace SeamPlacerImpl { -// In case there are custom enforcers/blockers, the loop polygon shall always have -// sides smaller than this (so it isn't limited to original resolution). -static constexpr float MINIMAL_POLYGON_SIDE = scaled(0.2f); - -// When spAligned is active and there is a support enforcer, -// add this penalty to its center. -static constexpr float ENFORCER_CENTER_PENALTY = -10.f; - - - -// This function was introduced in 2016 to assign penalties to overhangs. -// LukasM thinks that it discriminated a bit too much, so especially external -// seams were than placed in funny places (non-overhangs were preferred too much). -// He implemented his own version (below) which applies fixed penalty for really big overlaps. -// static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) -// { -// // The extrudate is not fully supported by the lower layer. Fit a polynomial penalty curve. -// // Solved by sympy package: -// /* -// from sympy import * -// (x,a,b,c,d,r,z)=symbols('x a b c d r z') -// p = a + b*x + c*x*x + d*x*x*x -// p2 = p.subs(solve([p.subs(x, -r), p.diff(x).subs(x, -r), p.diff(x,x).subs(x, -r), p.subs(x, 0)-z], [a, b, c, d])) -// from sympy.plotting import plot -// plot(p2.subs(r,0.2).subs(z,1.), (x, -1, 3), adaptive=False, nb_of_points=400) -// */ -// if (overlap_distance < - nozzle_r) { -// // The extrudate is fully supported by the lower layer. This is the ideal case, therefore zero penalty. -// return 0.f; -// } else { -// float x = overlap_distance / nozzle_r; -// float x2 = x * x; -// float x3 = x2 * x; -// return weight_zero * (1.f + 3.f * x + 3.f * x2 + x3); -// } -// } -static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) -{ - return overlap_distance > nozzle_r ? weight_zero : 0.f; +template int sgn(T val) { + return int(T(0) < val) - int(val < T(0)); } - - -// Return a value in <0, 1> of a cubic B-spline kernel centered around zero. -// The B-spline is re-scaled so it has value 1 at zero. -static inline float bspline_kernel(float x) -{ - x = std::abs(x); - if (x < 1.f) { - return 1.f - (3.f / 2.f) * x * x + (3.f / 4.f) * x * x * x; - } - else if (x < 2.f) { - x -= 1.f; - float x2 = x * x; - float x3 = x2 * x; - return (1.f / 4.f) - (3.f / 4.f) * x + (3.f / 4.f) * x2 - (1.f / 4.f) * x3; - } - else - return 0; +// base function: ((e^(((1)/(x^(2)+1)))-1)/(e-1)) +// checkout e.g. here: https://www.geogebra.org/calculator +float gauss(float value, float mean_x_coord, float mean_value, float falloff_speed) { + float shifted = value - mean_x_coord; + float denominator = falloff_speed * shifted * shifted + 1.0f; + float exponent = 1.0f / denominator; + return mean_value * (std::exp(exponent) - 1.0f) / (std::exp(1.0f) - 1.0f); } +/// Coordinate frame +class Frame { +public: + Frame() { + mX = Vec3f(1, 0, 0); + mY = Vec3f(0, 1, 0); + mZ = Vec3f(0, 0, 1); + } + Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) : + mX(x), mY(y), mZ(z) { + } -static Points::const_iterator project_point_to_polygon_and_insert(Polygon &polygon, const Point &pt, double eps) -{ - assert(polygon.points.size() >= 2); - if (polygon.points.size() <= 1) - if (polygon.points.size() == 1) - return polygon.points.begin(); + void set_from_z(const Vec3f &z) { + mZ = z.normalized(); + Vec3f tmpZ = mZ; + Vec3f tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3f(0, 1, 0) : Vec3f(1, 0, 0); + mY = (tmpZ.cross(tmpX)).normalized(); + mX = mY.cross(tmpZ); + } - Point pt_min; - double d_min = std::numeric_limits::max(); - size_t i_min = size_t(-1); + Vec3f to_world(const Vec3f &a) const { + return a.x() * mX + a.y() * mY + a.z() * mZ; + } - for (size_t i = 0; i < polygon.points.size(); ++ i) { - size_t j = i + 1; - if (j == polygon.points.size()) - j = 0; - const Point &p1 = polygon.points[i]; - const Point &p2 = polygon.points[j]; - const Slic3r::Point v_seg = p2 - p1; - const Slic3r::Point v_pt = pt - p1; - const int64_t l2_seg = int64_t(v_seg(0)) * int64_t(v_seg(0)) + int64_t(v_seg(1)) * int64_t(v_seg(1)); - int64_t t_pt = int64_t(v_seg(0)) * int64_t(v_pt(0)) + int64_t(v_seg(1)) * int64_t(v_pt(1)); - if (t_pt < 0) { - // Closest to p1. - double dabs = sqrt(int64_t(v_pt(0)) * int64_t(v_pt(0)) + int64_t(v_pt(1)) * int64_t(v_pt(1))); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - pt_min = p1; - } - } - else if (t_pt > l2_seg) { - // Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the next step. - continue; - } else { - // Closest to the segment. - assert(t_pt >= 0 && t_pt <= l2_seg); - int64_t d_seg = int64_t(v_seg(1)) * int64_t(v_pt(0)) - int64_t(v_seg(0)) * int64_t(v_pt(1)); - double d = double(d_seg) / sqrt(double(l2_seg)); - double dabs = std::abs(d); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - // Evaluate the foot point. - pt_min = p1; - double linv = double(d_seg) / double(l2_seg); - pt_min(0) = pt(0) - coord_t(floor(double(v_seg(1)) * linv + 0.5)); - pt_min(1) = pt(1) + coord_t(floor(double(v_seg(0)) * linv + 0.5)); - assert(Line(p1, p2).distance_to(pt_min) < scale_(1e-5)); - } + Vec3f to_local(const Vec3f &a) const { + return Vec3f(mX.dot(a), mY.dot(a), mZ.dot(a)); + } + + const Vec3f& binormal() const { + return mX; + } + + const Vec3f& tangent() const { + return mY; + } + + const Vec3f& normal() const { + return mZ; + } + +private: + Vec3f mX, mY, mZ; +}; + +Vec3f sample_sphere_uniform(const Vec2f &samples) { + float term1 = 2.0f * float(PI) * samples.x(); + float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); + return {cos(term1) * term2, sin(term1) * term2, + 1.0f - 2.0f * samples.y()}; +} + +Vec3f sample_hemisphere_uniform(const Vec2f &samples) { + float term1 = 2.0f * float(PI) * samples.x(); + float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); + return {cos(term1) * term2, sin(term1) * term2, + abs(1.0f - 2.0f * samples.y())}; +} + +Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) { + float term1 = 2.f * float(PI) * samples.x(); + float term2 = pow(samples.y(), 1.f / (power + 1.f)); + float term3 = sqrt(1.f - term2 * term2); + + return Vec3f(cos(term1) * term3, sin(term1) * term3, term2); +} + +std::vector raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree, + const indexed_triangle_set &triangles, size_t negative_volumes_start_index) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start"; + + //prepare uniform samples of a hemisphere + float step_size = 1.0f / SeamPlacer::sqr_rays_per_triangle; + std::vector precomputed_sample_directions( + SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); + for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_triangle; ++x_idx) { + float sample_x = x_idx * step_size + step_size / 2.0; + for (size_t y_idx = 0; y_idx < SeamPlacer::sqr_rays_per_triangle; ++y_idx) { + size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_triangle + y_idx; + float sample_y = y_idx * step_size + step_size / 2.0; + precomputed_sample_directions[dir_index] = sample_hemisphere_uniform( { sample_x, sample_y }); } } - assert(i_min != size_t(-1)); - if ((pt_min - polygon.points[i_min]).cast().norm() > eps) { - // Insert a new point on the segment i_min, i_min+1. - return polygon.points.insert(polygon.points.begin() + (i_min + 1), pt_min); - } - return polygon.points.begin() + i_min; -} + bool model_contains_negative_parts = negative_volumes_start_index < triangles.indices.size(); + std::vector result(triangles.indices.size()); + tbb::parallel_for(tbb::blocked_range(0, result.size()), + [&triangles, &precomputed_sample_directions, model_contains_negative_parts, negative_volumes_start_index, + &raycasting_tree, &result](tbb::blocked_range r) { + // Maintaining hits memory outside of the loop, so it does not have to be reallocated for each query. + std::vector hits; + for (size_t face_index = r.begin(); face_index < r.end(); ++face_index) { + FaceVisibilityInfo &dest = result[face_index]; + dest.visibility = 1.0f; + constexpr float decrease = 1.0f + / (SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); + Vec3i face = triangles.indices[face_index]; + Vec3f A = triangles.vertices[face.x()]; + Vec3f B = triangles.vertices[face.y()]; + Vec3f C = triangles.vertices[face.z()]; + Vec3f center = (A + B + C) / 3.0f; + Vec3f normal = ((B - A).cross(C - B)).normalized(); + // apply the local direction via Frame struct - the local_dir is with respect to +Z being forward + Frame f; + f.set_from_z(normal); -static std::vector polygon_angles_at_vertices(const Polygon &polygon, const std::vector &lengths, float min_arm_length) -{ - assert(polygon.points.size() + 1 == lengths.size()); - if (min_arm_length > 0.25f * lengths.back()) - min_arm_length = 0.25f * lengths.back(); - - // Find the initial prev / next point span. - size_t idx_prev = polygon.points.size(); - size_t idx_curr = 0; - size_t idx_next = 1; - while (idx_prev > idx_curr && lengths.back() - lengths[idx_prev] < min_arm_length) - -- idx_prev; - while (idx_next < idx_prev && lengths[idx_next] < min_arm_length) - ++ idx_next; - - std::vector angles(polygon.points.size(), 0.f); - for (; idx_curr < polygon.points.size(); ++ idx_curr) { - // Move idx_prev up until the distance between idx_prev and idx_curr is lower than min_arm_length. - if (idx_prev >= idx_curr) { - while (idx_prev < polygon.points.size() && lengths.back() - lengths[idx_prev] + lengths[idx_curr] > min_arm_length) - ++ idx_prev; - if (idx_prev == polygon.points.size()) - idx_prev = 0; - } - while (idx_prev < idx_curr && lengths[idx_curr] - lengths[idx_prev] > min_arm_length) - ++ idx_prev; - // Move idx_prev one step back. - if (idx_prev == 0) - idx_prev = polygon.points.size() - 1; - else - -- idx_prev; - // Move idx_next up until the distance between idx_curr and idx_next is greater than min_arm_length. - if (idx_curr <= idx_next) { - while (idx_next < polygon.points.size() && lengths[idx_next] - lengths[idx_curr] < min_arm_length) - ++ idx_next; - if (idx_next == polygon.points.size()) - idx_next = 0; - } - while (idx_next < idx_curr && lengths.back() - lengths[idx_curr] + lengths[idx_next] < min_arm_length) - ++ idx_next; - // Calculate angle between idx_prev, idx_curr, idx_next. - const Point &p0 = polygon.points[idx_prev]; - const Point &p1 = polygon.points[idx_curr]; - const Point &p2 = polygon.points[idx_next]; - const Point v1 = p1 - p0; - const Point v2 = p2 - p1; - int64_t dot = int64_t(v1(0))*int64_t(v2(0)) + int64_t(v1(1))*int64_t(v2(1)); - int64_t cross = int64_t(v1(0))*int64_t(v2(1)) - int64_t(v1(1))*int64_t(v2(0)); - float angle = float(atan2(double(cross), double(dot))); - angles[idx_curr] = angle; - } - - return angles; -} - - - -void SeamPlacer::init(const Print& print) -{ - m_enforcers.clear(); - m_blockers.clear(); - m_seam_history.clear(); - m_po_list.clear(); - - const std::vector& nozzle_dmrs = print.config().nozzle_diameter.values; - float max_nozzle_dmr = *std::max_element(nozzle_dmrs.begin(), nozzle_dmrs.end()); - - - std::vector temp_enf; - std::vector temp_blk; - std::vector temp_polygons; - - for (const PrintObject* po : print.objects()) { - - auto merge_and_offset = [po, &temp_polygons, max_nozzle_dmr](EnforcerBlockerType type, std::vector& out) { - // Offset the triangles out slightly. - auto offset_out = [](Polygon& input, float offset) -> ExPolygons { - ClipperLib::Paths out(1); - std::vector deltas(input.points.size(), offset); - input.make_counter_clockwise(); - out.front() = mittered_offset_path_scaled(input.points, deltas, 3.); - return ClipperPaths_to_Slic3rExPolygons(out, true); // perform union - }; - - - temp_polygons.clear(); - po->project_and_append_custom_facets(true, type, temp_polygons); - out.clear(); - out.reserve(temp_polygons.size()); - float offset = scale_(max_nozzle_dmr + po->config().elefant_foot_compensation); - for (Polygons &src : temp_polygons) { - out.emplace_back(ExPolygons()); - for (Polygon& plg : src) { - ExPolygons offset_explg = offset_out(plg, offset); - if (! offset_explg.empty()) - out.back().emplace_back(std::move(offset_explg.front())); - } - - offset = scale_(max_nozzle_dmr); - } - }; - merge_and_offset(EnforcerBlockerType::BLOCKER, temp_blk); - merge_and_offset(EnforcerBlockerType::ENFORCER, temp_enf); - - // Remember this PrintObject and initialize a store of enforcers and blockers for it. - m_po_list.push_back(po); - size_t po_idx = m_po_list.size() - 1; - m_enforcers.emplace_back(std::vector(temp_enf.size())); - m_blockers.emplace_back(std::vector(temp_blk.size())); - - // A helper class to store data to build the AABB tree from. - class CustomTriangleRef { - public: - CustomTriangleRef(size_t idx, - Point&& centroid, - BoundingBox&& bb) - : m_idx{idx}, m_centroid{centroid}, - m_bbox{AlignedBoxType(bb.min, bb.max)} - {} - size_t idx() const { return m_idx; } - const Point& centroid() const { return m_centroid; } - const TreeType::BoundingBox& bbox() const { return m_bbox; } - - private: - size_t m_idx; - Point m_centroid; - AlignedBoxType m_bbox; - }; - - // A lambda to extract the ExPolygons and save them into the member AABB tree. - // Will be called for enforcers and blockers separately. - auto add_custom = [](std::vector& src, std::vector& dest) { - // Go layer by layer, and append all the ExPolygons into the AABB tree. - size_t layer_idx = 0; - for (ExPolygons& expolys_on_layer : src) { - CustomTrianglesPerLayer& layer_data = dest[layer_idx]; - std::vector triangles_data; - layer_data.polys.reserve(expolys_on_layer.size()); - triangles_data.reserve(expolys_on_layer.size()); - - for (ExPolygon& expoly : expolys_on_layer) { - if (expoly.empty()) - continue; - layer_data.polys.emplace_back(std::move(expoly)); - triangles_data.emplace_back(layer_data.polys.size() - 1, - layer_data.polys.back().centroid(), - layer_data.polys.back().bounding_box()); - } - // All polygons are saved, build the AABB tree for them. - layer_data.tree.build(std::move(triangles_data)); - ++layer_idx; - } - }; - - add_custom(temp_enf, m_enforcers.at(po_idx)); - add_custom(temp_blk, m_blockers.at(po_idx)); - } -} - - - -void SeamPlacer::plan_perimeters(const std::vector perimeters, - const Layer& layer, SeamPosition seam_position, - Point last_pos, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid) -{ - // When printing the perimeters, we want the seams on external and internal perimeters to match. - // We have a list of perimeters in the order to be printed. Each internal perimeter must inherit - // the seam from the previous external perimeter. - - m_plan.clear(); - m_plan_idx = 0; - - if (perimeters.empty() || ! po) - return; - - m_plan.resize(perimeters.size()); - - for (int i = 0; i < int(perimeters.size()); ++i) { - if (perimeters[i]->role() == erExternalPerimeter && perimeters[i]->is_loop()) { - last_pos = this->calculate_seam( - layer, seam_position, *dynamic_cast(perimeters[i]), nozzle_dmr, - po, lower_layer_edge_grid, last_pos, false); - m_plan[i].external = true; - } - m_plan[i].seam_position = seam_position; - m_plan[i].layer = &layer; - m_plan[i].po = po; - m_plan[i].pt = last_pos; - } -} - - -void SeamPlacer::place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter, - const EdgeGrid::Grid* lower_layer_edge_grid) -{ - // const double seam_offset = nozzle_diameter; - - Point seam = last_pos; - if (! m_plan.empty() && m_plan_idx < m_plan.size()) { - if (m_plan[m_plan_idx].external) { - seam = m_plan[m_plan_idx].pt; - // One more heuristics: if the seam is too far from current nozzle position, - // try to place it again. This can happen in cases where the external perimeter - // does not belong to the preceding ones and they are ordered so they end up - // far from each other. - if ((seam.cast() - last_pos.cast()).squaredNorm() > std::pow(scale_(5.*nozzle_diameter), 2.)) - seam = this->calculate_seam(*m_plan[m_plan_idx].layer, m_plan[m_plan_idx].seam_position, loop, nozzle_diameter, - m_plan[m_plan_idx].po, lower_layer_edge_grid, last_pos, false); - - if (m_plan[m_plan_idx].seam_position == spAligned) - m_seam_history.add_seam(m_plan[m_plan_idx].po, m_plan[m_plan_idx].pt, loop.polygon().bounding_box()); - } - else { - if (!external_first) { - // Internal perimeter printed before the external. - // First get list of external seams. - std::vector ext_seams; - size_t external_cnt = 0; - for (size_t i = 0; i < m_plan.size(); ++i) { - if (m_plan[i].external) { - ext_seams.emplace_back(i); - ++external_cnt; - } - } - - if (!ext_seams.empty()) { - // First find the line segment closest to an external seam: - //int path_idx = 0; - //int line_idx = 0; - size_t ext_seam_idx = size_t(-1); - double min_dist_sqr = std::numeric_limits::max(); - std::vector lines_vect; - for (int i = 0; i < int(loop.paths.size()); ++i) { - lines_vect.emplace_back(loop.paths[i].polyline.lines()); - const Lines& lines = lines_vect.back(); - for (int j = 0; j < int(lines.size()); ++j) { - for (size_t k : ext_seams) { - double d_sqr = lines[j].distance_to_squared(m_plan[k].pt); - if (d_sqr < min_dist_sqr) { - //path_idx = i; - //line_idx = j; - ext_seam_idx = k; - min_dist_sqr = d_sqr; + for (const auto &dir : precomputed_sample_directions) { + Vec3f final_ray_dir = (f.to_world(dir)); + if (!model_contains_negative_parts) { + igl::Hit hitpoint; + // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and + // direction. + Vec3d final_ray_dir_d = final_ray_dir.cast(); + Vec3d ray_origin_d = (center + normal * 0.1).cast(); // start above surface. + bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices, + triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint); + if (hit) { + dest.visibility -= decrease; + } + } else { //TODO improve logic for order based boolean operations - consider order of volumes + Vec3d ray_origin_d = (center + normal * 0.1).cast(); // start above surface. + if (face_index >= negative_volumes_start_index) { // if casting from negative volume face, invert direction, change start pos + final_ray_dir = -1.0 * final_ray_dir; + ray_origin_d = (center - normal * 0.1).cast(); + } + Vec3d final_ray_dir_d = final_ray_dir.cast(); + bool some_hit = AABBTreeIndirect::intersect_ray_all_hits(triangles.vertices, + triangles.indices, raycasting_tree, + ray_origin_d, final_ray_dir_d, hits); + if (some_hit) { + int in_negative = 0; + int in_positive = 0; + // NOTE: iterating in reverse, from the last hit for one simple reason: We know the state of the ray at that point; + // It cannot be inside model, and it cannot be inside negative volume + for (int hit_index = int(hits.size()) - 1; hit_index >= 0; --hit_index) { + Vec3f face_normal = its_face_normal(triangles, hits[hit_index].id); + if (hits[hit_index].id >= int(negative_volumes_start_index)) { //negative volume hit + in_negative += sgn(face_normal.dot(final_ray_dir)); // if volume face aligns with ray dir, we are leaving negative space + // which in reverse hit analysis means, that we are entering negative space :) and vice versa + } else { + in_positive += sgn(face_normal.dot(final_ray_dir)); + } + if (in_positive > 0 && in_negative <= 0) { + dest.visibility -= decrease; + break; + } } } } } + } + }); - // Only accept seam that is reasonably close. - if (ext_seam_idx != size_t(-1)) { - // How many nozzle diameters is considered "close"? - const double nozzle_d_limit = 2. * (1. + m_plan.size() / external_cnt); - const double limit_dist_sqr = double(scale_(scale_((unscale(m_plan[ext_seam_idx].pt) - unscale(m_plan[m_plan_idx].pt)).squaredNorm() * std::pow(nozzle_d_limit * nozzle_diameter, 2.)))); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end"; - if (min_dist_sqr < limit_dist_sqr) { - // Now find a projection of the external seam - //const Lines& lines = lines_vect[path_idx]; - //Point closest = m_plan[ext_seam_idx].pt.projection_onto(lines[line_idx]); + return result; +} - // This code does staggering of internal perimeters, turned off for now. - // - // double dist = (closest.cast() - lines[line_idx].b.cast()).norm(); - // - // // And walk along the perimeter until we make enough space for - // // seams of all perimeters beforethe external one. - // double offset = (ext_seam_idx - m_plan_idx) * scale_(seam_offset); - // double last_offset = offset; - // offset -= dist; - // const Point* a = &closest; - // const Point* b = &lines[line_idx].b; - // while (++line_idx < int(lines.size()) && offset > 0.) { - // last_offset = offset; - // offset -= lines[line_idx].length(); - // a = &lines[line_idx].a; - // b = &lines[line_idx].b; - // } - // - // // We have walked far enough, too far maybe. Interpolate on the - // // last segment to find the end precisely. - // offset = std::min(0., offset); // In case that offset is still positive (we may have "wrapped around") - // double ratio = last_offset / (last_offset - offset); - // seam = (a->cast() + ((b->cast() - a->cast()) * ratio)).cast(); - seam = m_plan[ext_seam_idx].pt; - } +std::vector calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector &lengths, + float min_arm_length) { + std::vector result(polygon.size()); + + if (polygon.size() == 1) { + result[0] = 0.0f; + } + + size_t idx_prev = 0; + size_t idx_curr = 0; + size_t idx_next = 0; + + float distance_to_prev = 0; + float distance_to_next = 0; + + //push idx_prev far enough back as initialization + while (distance_to_prev < min_arm_length) { + idx_prev = Slic3r::prev_idx_modulo(idx_prev, polygon.size()); + distance_to_prev += lengths[idx_prev]; + } + + for (size_t _i = 0; _i < polygon.size(); ++_i) { + // pull idx_prev to current as much as possible, while respecting the min_arm_length + while (distance_to_prev - lengths[idx_prev] > min_arm_length) { + distance_to_prev -= lengths[idx_prev]; + idx_prev = Slic3r::next_idx_modulo(idx_prev, polygon.size()); + } + + //push idx_next forward as far as needed + while (distance_to_next < min_arm_length) { + distance_to_next += lengths[idx_next]; + idx_next = Slic3r::next_idx_modulo(idx_next, polygon.size()); + } + + // Calculate angle between idx_prev, idx_curr, idx_next. + const Point &p0 = polygon.points[idx_prev]; + const Point &p1 = polygon.points[idx_curr]; + const Point &p2 = polygon.points[idx_next]; + result[idx_curr] = float(angle(p1 - p0, p2 - p1)); + + // increase idx_curr by one + float curr_distance = lengths[idx_curr]; + idx_curr++; + distance_to_prev += curr_distance; + distance_to_next -= curr_distance; + } + + return result; +} + +// structure to store global information about the model - occlusion hits, enforcers, blockers +struct GlobalModelInfo { + indexed_triangle_set model; + AABBTreeIndirect::Tree<3, float> model_tree; + std::vector visiblity_info; + indexed_triangle_set enforcers; + indexed_triangle_set blockers; + AABBTreeIndirect::Tree<3, float> enforcers_tree; + AABBTreeIndirect::Tree<3, float> blockers_tree; + + bool is_enforced(const Vec3f &position, float radius) const { + if (enforcers.empty()) { + return false; + } + float radius_sqr = radius * radius; + return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices, + enforcers_tree, position, radius_sqr); + } + + bool is_blocked(const Vec3f &position, float radius) const { + if (blockers.empty()) { + return false; + } + float radius_sqr = radius * radius; + return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices, + blockers_tree, position, radius_sqr); + } + + float calculate_point_visibility(const Vec3f &position) const { + size_t hit_idx; + Vec3f hit_point; + if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree, + position, hit_idx, hit_point) >= 0) { + return visiblity_info[hit_idx].visibility; + } else { + return 0.0f; + } + + } + +#ifdef DEBUG_FILES + void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const { + indexed_triangle_set divided_mesh = obj_mesh; + Slic3r::CNumericLocalesSetter locales_setter; + + FILE *fp = boost::nowide::fopen(file_name, "w"); + if (fp == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << file_name << " for writing"; + return; + } + + for (size_t i = 0; i < divided_mesh.vertices.size(); ++i) { + float visibility = calculate_point_visibility(divided_mesh.vertices[i]); + Vec3f color = value_to_rgbf(0.0f, 1.0f, + visibility); + fprintf(fp, "v %f %f %f %f %f %f\n", + divided_mesh.vertices[i](0), divided_mesh.vertices[i](1), divided_mesh.vertices[i](2), + color(0), color(1), color(2) + ); + } + for (size_t i = 0; i < divided_mesh.indices.size(); ++i) + fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1, + divided_mesh.indices[i][2] + 1); + fclose(fp); + + } +#endif +} +; + +//Extract perimeter polygons of the given layer +Polygons extract_perimeter_polygons(const Layer *layer, const SeamPosition configured_seam_preference, + std::vector &corresponding_regions_out) { + Polygons polygons; + for (const LayerRegion *layer_region : layer->regions()) { + for (const ExtrusionEntity *ex_entity : layer_region->perimeters.entities) { + if (ex_entity->is_collection()) { //collection of inner, outer, and overhang perimeters + for (const ExtrusionEntity *perimeter : static_cast(ex_entity)->entities) { + if (perimeter->role() == ExtrusionRole::erExternalPerimeter + || (perimeter->role() == ExtrusionRole::erPerimeter + && configured_seam_preference == spRandom)) { //for random seam alignment, extract all perimeters + Points p; + perimeter->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); } } + if (polygons.empty()) { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); + } + } else { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); } - else { - // We should have a candidate ready from before. If not, use last_pos. - if (m_plan_idx > 0 && m_plan[m_plan_idx - 1].precalculated) - seam = m_plan[m_plan_idx - 1].pt; - } - - // seam now contains a hot candidate for internal seam. Use it unless there is a sharp corner nearby. - // We will call the normal seam planning function, pretending that we are currently at the candidate point - // and set to spNearest. If the ideal seam it finds is close to current candidate, use it. - // This is to prevent having seams very close to corners, just because of external seam position. - seam = calculate_seam(*m_plan[m_plan_idx].layer, spNearest, loop, nozzle_diameter, - m_plan[m_plan_idx].po, lower_layer_edge_grid, seam, true); } - m_plan[m_plan_idx].pt = seam; } - - // Split the loop at the point with a minium penalty. - if (!loop.split_at_vertex(seam)) - // The point is not in the original loop. Insert it. - loop.split_at(seam, true); - - if (external_first && m_plan_idx+1 1) { -// const ExtrusionPath& last = loop.paths.back(); -// auto it = last.polyline.points.crbegin() + 1; -// for (; it != last.polyline.points.crend(); ++it) { -// running_sqr += (it->cast() - (it - 1)->cast()).squaredNorm(); -// if (running_sqr > dist_sqr) -// break; -// running_sqr_last = running_sqr; -// } -// if (running_sqr <= dist_sqr) -// it = last.polyline.points.crend() - 1; -// // Now interpolate. -// double ratio = (std::sqrt(dist_sqr) - std::sqrt(running_sqr_last)) / (std::sqrt(running_sqr) - std::sqrt(running_sqr_last)); -// m_plan[m_plan_idx + 1].pt = ((it - 1)->cast() + (it->cast() - (it - 1)->cast()) * std::min(ratio, 1.)).cast(); -// m_plan[m_plan_idx + 1].precalculated = true; - m_plan[m_plan_idx + 1].pt = m_plan[m_plan_idx].pt; - m_plan[m_plan_idx + 1].precalculated = true; -// } + if (polygons.empty()) { // If there are no perimeter polygons for whatever reason (disabled perimeters .. ) insert dummy point + // it is easier than checking everywhere if the layer is not emtpy, no seam will be placed to this layer anyway + polygons.emplace_back(std::vector { Point { 0, 0 } }); + corresponding_regions_out.push_back(nullptr); } - ++m_plan_idx; + return polygons; } +// Insert SeamCandidates created from perimeter polygons in to the result vector. +// Compute its type (Enfrocer,Blocker), angle, and position +//each SeamCandidate also contains pointer to shared Perimeter structure representing the polygon +// if Custom Seam modifiers are present, oversamples the polygon if necessary to better fit user intentions +void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const LayerRegion *region, + const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) { + if (orig_polygon.size() == 0) { + return; + } -// Returns "best" seam for a given perimeter. -Point SeamPlacer::calculate_seam(const Layer& layer, const SeamPosition seam_position, - const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest) -{ - Polygon polygon = loop.polygon(); + Polygon polygon = orig_polygon; bool was_clockwise = polygon.make_counter_clockwise(); - const coord_t nozzle_r = coord_t(scale_(0.5 * nozzle_dmr) + 0.5); - size_t po_idx = std::find(m_po_list.begin(), m_po_list.end(), po) - m_po_list.begin(); - - // Find current layer in respective PrintObject. Cache the result so the - // lookup is only done once per layer, not for each loop. - const Layer* layer_po = nullptr; - if (po == m_last_po && layer.print_z == m_last_print_z) - layer_po = m_last_layer_po; - else { - layer_po = po ? po->get_layer_at_printz(layer.print_z) : nullptr; - m_last_po = po; - m_last_print_z = layer.print_z; - m_last_layer_po = layer_po; + std::vector lengths { }; + for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) { + lengths.push_back(std::max((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm(), 0.01)); } - if (! layer_po) - return last_pos; + lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.01)); - // Index of this layer in the respective PrintObject. - size_t layer_idx = layer_po->id() - po->layers().front()->id(); // raft layers + std::vector local_angles = calculate_polygon_angles_at_vertices(polygon, lengths, + SeamPlacer::polygon_local_angles_arm_distance); - assert(layer_idx < po->layer_count()); + result.perimeters.push_back( { }); + Perimeter &perimeter = result.perimeters.back(); - const bool custom_seam = loop.role() == erExternalPerimeter && this->is_custom_seam_on_layer(layer_idx, po_idx); - - if (custom_seam) { - // Seam enf/blockers can begin and end in between the original vertices. - // Let add extra points in between and update the leghths. - polygon.densify(MINIMAL_POLYGON_SIDE); + std::queue orig_polygon_points { }; + for (size_t index = 0; index < polygon.size(); ++index) { + Vec2f unscaled_p = unscale(polygon[index]).cast(); + orig_polygon_points.emplace(unscaled_p.x(), unscaled_p.y(), z_coord); } + Vec3f first = orig_polygon_points.front(); + std::queue oversampled_points { }; + size_t orig_angle_index = 0; + perimeter.start_index = result.points.size(); + perimeter.flow_width = region != nullptr ? region->flow(FlowRole::frExternalPerimeter).width() : 0.0f; + bool some_point_enforced = false; + while (!orig_polygon_points.empty() || !oversampled_points.empty()) { + EnforcedBlockedSeamPoint type = EnforcedBlockedSeamPoint::Neutral; + Vec3f position; + float local_ccw_angle = 0; + bool orig_point = false; + if (!oversampled_points.empty()) { + position = oversampled_points.front(); + oversampled_points.pop(); + } else { + position = orig_polygon_points.front(); + orig_polygon_points.pop(); + local_ccw_angle = was_clockwise ? -local_angles[orig_angle_index] : local_angles[orig_angle_index]; + orig_angle_index++; + orig_point = true; + } - if (seam_position != spRandom) { - // Retrieve the last start position for this object. - float last_pos_weight = 1.f; + if (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { + type = EnforcedBlockedSeamPoint::Enforced; + some_point_enforced = true; + } - if (seam_position == spAligned) { - // Seam is aligned to the seam at the preceding layer. - if (po != nullptr) { - std::optional pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, loop.polygon().bounding_box()); - if (pos.has_value()) { - last_pos = *pos; - last_pos_weight = is_custom_enforcer_on_layer(layer_idx, po_idx) ? 0.f : 1.f; + if (global_model_info.is_blocked(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { + type = EnforcedBlockedSeamPoint::Blocked; + } + + if (orig_point) { + Vec3f pos_of_next = orig_polygon_points.empty() ? first : orig_polygon_points.front(); + float distance_to_next = (position - pos_of_next).norm(); + if (global_model_info.is_enforced(position, distance_to_next)) { + Vec3f vec_to_next = (pos_of_next - position).normalized(); + float step_size = SeamPlacer::enforcer_oversampling_distance; + float step = step_size; + while (step < distance_to_next) { + oversampled_points.push(position + vec_to_next * step); + step += step_size; } } } - else if (seam_position == spRear) { - // Object is centered around (0,0) in its current coordinate system. - last_pos.x() = 0; - last_pos.y() = coord_t(3. * po->bounding_box().radius()); - last_pos_weight = 5.f; - } if (seam_position == spNearest) { - // last_pos already contains current nozzle position - } - // Insert a projection of last_pos into the polygon. - size_t last_pos_proj_idx; - { - auto it = project_point_to_polygon_and_insert(polygon, last_pos, 0.1 * nozzle_r); - last_pos_proj_idx = it - polygon.points.begin(); - } - - // Parametrize the polygon by its length. - std::vector lengths = polygon.parameter_by_length(); - - // For each polygon point, store a penalty. - // First calculate the angles, store them as penalties. The angles are caluculated over a minimum arm length of nozzle_r. - std::vector penalties = polygon_angles_at_vertices(polygon, lengths, float(nozzle_r)); - // No penalty for reflex points, slight penalty for convex points, high penalty for flat surfaces. - const float penaltyConvexVertex = 1.f; - const float penaltyFlatSurface = 3.f; - const float penaltyOverhangHalf = 10.f; - // Penalty for visible seams. - for (size_t i = 0; i < polygon.points.size(); ++ i) { - float ccwAngle = penalties[i]; - if (was_clockwise) - ccwAngle = - ccwAngle; - float penalty = 0; - if (ccwAngle <- float(0.6 * PI)) - // Sharp reflex vertex. We love that, it hides the seam perfectly. - penalty = 0.f; - else if (ccwAngle > float(0.6 * PI)) - // Seams on sharp convex vertices are more visible than on reflex vertices. - penalty = penaltyConvexVertex; - else if (ccwAngle < 0.f) { - // Interpolate penalty between maximum and zero. - penalty = penaltyFlatSurface * bspline_kernel(ccwAngle * float(PI * 2. / 3.)); - } else { - assert(ccwAngle >= 0.f); - // Interpolate penalty between maximum and the penalty for a convex vertex. - penalty = penaltyConvexVertex + (penaltyFlatSurface - penaltyConvexVertex) * bspline_kernel(ccwAngle * float(PI * 2. / 3.)); - } - // Give a negative penalty for points close to the last point or the prefered seam location. - float dist_to_last_pos_proj = (i < last_pos_proj_idx) ? - std::min(lengths[last_pos_proj_idx] - lengths[i], lengths.back() - lengths[last_pos_proj_idx] + lengths[i]) : - std::min(lengths[i] - lengths[last_pos_proj_idx], lengths.back() - lengths[i] + lengths[last_pos_proj_idx]); - float dist_max = 0.1f * lengths.back(); // 5.f * nozzle_dmr - penalty -= last_pos_weight * bspline_kernel(dist_to_last_pos_proj / dist_max); - penalties[i] = std::max(0.f, penalty); - if (prefer_nearest) { - // This hack limits the search around the nearest position projection. - penalties[i] += dist_to_last_pos_proj > 6.f * nozzle_r ? 100.f : 0.f; - } - } - - // Penalty for overhangs. - if (lower_layer_edge_grid) { - // Use the edge grid distance field structure over the lower layer to calculate overhangs. - coord_t nozzle_r = coord_t(std::floor(scale_(0.5 * nozzle_dmr) + 0.5)); - coord_t search_r = coord_t(std::floor(scale_(0.8 * nozzle_dmr) + 0.5)); - for (size_t i = 0; i < polygon.points.size(); ++ i) { - const Point &p = polygon.points[i]; - coordf_t dist; - // Signed distance is positive outside the object, negative inside the object. - // The point is considered at an overhang, if it is more than nozzle radius - // outside of the lower layer contour. - [[maybe_unused]] bool found = lower_layer_edge_grid->signed_distance(p, search_r, dist); - // If the approximate Signed Distance Field was initialized over lower_layer_edge_grid, - // then the signed distnace shall always be known. - assert(found); - penalties[i] += extrudate_overlap_penalty(float(nozzle_r), penaltyOverhangHalf, float(dist)); - } - } - - // Custom seam. Huge (negative) constant penalty is applied inside - // blockers (enforcers) to rule out points that should not win. - if (custom_seam) - this->apply_custom_seam(polygon, po_idx, penalties, lengths, layer_idx, seam_position); - - // Find a point with a minimum penalty. - size_t idx_min = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); - - if (seam_position != spAligned || ! is_custom_enforcer_on_layer(layer_idx, po_idx)) { - // Very likely the weight of idx_min is very close to the weight of last_pos_proj_idx. - // In that case use last_pos_proj_idx instead. - float penalty_aligned = penalties[last_pos_proj_idx]; - float penalty_min = penalties[idx_min]; - float penalty_diff_abs = std::abs(penalty_min - penalty_aligned); - float penalty_max = std::max(std::abs(penalty_min), std::abs(penalty_aligned)); - float penalty_diff_rel = (penalty_max == 0.f) ? 0.f : penalty_diff_abs / penalty_max; - // printf("Align seams, penalty aligned: %f, min: %f, diff abs: %f, diff rel: %f\n", penalty_aligned, penalty_min, penalty_diff_abs, penalty_diff_rel); - if (std::abs(penalty_diff_rel) < 0.05) { - // Penalty of the aligned point is very close to the minimum penalty. - // Align the seams as accurately as possible. - idx_min = last_pos_proj_idx; - } - } - - - // Export the contour into a SVG file. - #if 0 - { - static int iRun = 0; - SVG svg(debug_out_path("GCode_extrude_loop-%d.svg", iRun ++)); - if (m_layer->lower_layer != NULL) - svg.draw(m_layer->lower_layer->slices); - for (size_t i = 0; i < loop.paths.size(); ++ i) - svg.draw(loop.paths[i].as_polyline(), "red"); - Polylines polylines; - for (size_t i = 0; i < loop.paths.size(); ++ i) - polylines.push_back(loop.paths[i].as_polyline()); - Slic3r::Polygons polygons; - coordf_t nozzle_dmr = EXTRUDER_CONFIG(nozzle_diameter); - coord_t delta = scale_(0.5*nozzle_dmr); - Slic3r::offset(polylines, &polygons, delta); -// for (size_t i = 0; i < polygons.size(); ++ i) svg.draw((Polyline)polygons[i], "blue"); - svg.draw(last_pos, "green", 3); - svg.draw(polygon.points[idx_min], "yellow", 3); - svg.Close(); - } - #endif - return polygon.points[idx_min]; - - } else - return this->get_random_seam(layer_idx, polygon, po_idx); -} - - -Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_idx, - bool* saw_custom) const -{ - // Parametrize the polygon by its length. - const std::vector lengths = polygon.parameter_by_length(); - - // Which of the points are inside enforcers/blockers? - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_idx, polygon, po_idx, enforcers_idxs, blockers_idxs); - - bool has_enforcers = ! enforcers_idxs.empty(); - bool has_blockers = ! blockers_idxs.empty(); - if (saw_custom) - *saw_custom = has_enforcers || has_blockers; - - assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end())); - assert(std::is_sorted(blockers_idxs.begin(), blockers_idxs.end())); - std::vector edges; - - // Lambda to calculate lengths of all edges of interest. Last parameter - // decides whether to measure edges inside or outside idxs. - // Negative number = not an edge of interest. - auto get_valid_length = [&lengths](const std::vector& idxs, - std::vector& edges, - bool measure_inside_edges) -> float - { - // First mark edges we are interested in by assigning a positive number. - edges.assign(lengths.size()-1, measure_inside_edges ? -1.f : 1.f); - for (size_t i=0; i the edge between them is the enforcer/blocker. - bool inside_edge = ((i != idxs.size()-1 && idxs[i+1] == this_pt_idx + 1) - || (i == idxs.size()-1 && idxs.back() == lengths.size()-2 && idxs[0] == 0)); - if (inside_edge) - edges[this_pt_idx] = measure_inside_edges ? 1.f : -1.f; - } - // Now measure them. - float running_total = 0.f; - for (size_t i=0; i 0.f) { - edges[i] = lengths[i+1] - lengths[i]; - running_total += edges[i]; - } - } - return running_total; - }; - - // Find all seam candidate edges and their lengths. - float valid_length = 0.f; - if (has_enforcers) - valid_length = get_valid_length(enforcers_idxs, edges, true); - - if (! has_enforcers || valid_length == 0.f) { - // Second condition covers case with isolated enf points. Given how the painted - // triangles are projected, this should not happen. Stay on the safe side though. - if (has_blockers) - valid_length = get_valid_length(blockers_idxs, edges, false); - if (valid_length == 0.f) // No blockers or everything blocked - use the whole polygon. - valid_length = lengths.back(); + result.points.emplace_back(position, perimeter, local_ccw_angle, type); } - assert(valid_length != 0.f); - // Now generate a random length and find the respective edge. - float rand_len = valid_length * (rand()/float(RAND_MAX)); - size_t pt_idx = 0; // Index of the edge where to put the seam. - if (valid_length == lengths.back()) { - // Whole polygon is used for placing the seam. - auto it = std::lower_bound(lengths.begin(), lengths.end(), rand_len); - pt_idx = it == lengths.begin() ? 0 : (it-lengths.begin()-1); // this takes care of a corner case where rand() returns 0 - } else { - float running = 0.f; - for (size_t i=0; i 0.f ? edges[i] : 0.f; - if (running >= rand_len) { - pt_idx = i; + + perimeter.end_index = result.points.size() - 1; + + // We will find first patch of enforced points (patch: continuous section of enforced points) and select the middle + // point, which will have priority during alignment + // If there are multiple enforced patches in the perimeter, others are ignored + if (some_point_enforced) { + size_t perimeter_size = perimeter.end_index - perimeter.start_index + 1; + const auto next_index = [&](size_t idx) { + return perimeter.start_index + Slic3r::next_idx_modulo(idx - perimeter.start_index, perimeter_size); + }; + + size_t first_enforced_idx = perimeter.start_index; + for (size_t _ = 0; _ < perimeter_size; ++_) { + if (result.points[first_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced && + result.points[next_index(first_enforced_idx)].type == EnforcedBlockedSeamPoint::Enforced) { break; } + first_enforced_idx = next_index(first_enforced_idx); } - } + first_enforced_idx = next_index(first_enforced_idx); - if (! has_enforcers && ! has_blockers) { - // The polygon may be too coarse, calculate the point exactly. - assert(valid_length == lengths.back()); - bool last_seg = pt_idx == polygon.points.size()-1; - size_t next_idx = last_seg ? 0 : pt_idx+1; - const Point& prev = polygon.points[pt_idx]; - const Point& next = polygon.points[next_idx]; - assert(next_idx == 0 || pt_idx+1 == next_idx); - coordf_t diff_x = next.x() - prev.x(); - coordf_t diff_y = next.y() - prev.y(); - coordf_t dist = lengths[last_seg ? pt_idx+1 : next_idx] - lengths[pt_idx]; - return Point(prev.x() + (rand_len - lengths[pt_idx]) * (diff_x/dist), - prev.y() + (rand_len - lengths[pt_idx]) * (diff_y/dist)); - - } else { - // The polygon should be dense enough. - return polygon.points[pt_idx]; - } -} - - - - - - - - -void SeamPlacer::get_enforcers_and_blockers(size_t layer_id, - const Polygon& polygon, - size_t po_idx, - std::vector& enforcers_idxs, - std::vector& blockers_idxs) const -{ - enforcers_idxs.clear(); - blockers_idxs.clear(); - - auto is_inside = [](const Point& pt, - const CustomTrianglesPerLayer& custom_data) -> bool { - assert(! custom_data.polys.empty()); - // Now ask the AABB tree which polygons we should check and check them. - std::vector candidates; - AABBTreeIndirect::get_candidate_idxs(custom_data.tree, pt, candidates); - if (! candidates.empty()) - for (size_t idx : candidates) - if (custom_data.polys[idx].contains(pt)) - return true; - return false; - }; - - if (! m_enforcers[po_idx].empty()) { - const CustomTrianglesPerLayer& enforcers = m_enforcers[po_idx][layer_id]; - if (! enforcers.polys.empty()) { - for (size_t i=0; i orig_large_angle_points_indices { }; + std::vector viable_points_indices { }; + size_t last_enforced_idx = first_enforced_idx; + for (size_t _ = 0; _ < perimeter_size; ++_) { + if (result.points[last_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced) { + break; } + viable_points_indices.push_back(last_enforced_idx); + if (abs(result.points[last_enforced_idx].local_ccw_angle) > 0.4 * PI) { + orig_large_angle_points_indices.push_back(last_enforced_idx); + } + last_enforced_idx = next_index(last_enforced_idx); } - } - - if (! m_blockers[po_idx].empty()) { - const CustomTrianglesPerLayer& blockers = m_blockers[po_idx][layer_id]; - if (! blockers.polys.empty()) { - for (size_t i=0; i find_enforcer_centers(const Polygon& polygon, - const std::vector& lengths, - const std::vector& enforcers_idxs) -{ - std::vector out; - assert(polygon.points.size()+1 == lengths.size()); - assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end())); - if (polygon.size() < 2 || enforcers_idxs.empty()) - return out; - - auto get_center_idx = [&lengths](size_t start_idx, size_t end_idx) -> size_t { - assert(end_idx >= start_idx); - if (start_idx == end_idx) - return start_idx; - float t_c = lengths[start_idx] + 0.5f * (lengths[end_idx] - lengths[start_idx]); - auto it = std::lower_bound(lengths.begin() + start_idx, lengths.begin() + end_idx, t_c); - int ret = it - lengths.begin(); - return ret; - }; - - int last_enforcer_start_idx = enforcers_idxs.front(); - bool first_pt_in_list = enforcers_idxs.front() != 0; - bool last_pt_in_list = enforcers_idxs.back() == polygon.points.size() - 1; - bool wrap_around = last_pt_in_list && first_pt_in_list; - - for (size_t i=0; i 0); + if (orig_large_angle_points_indices.empty()) { + size_t central_idx = viable_points_indices[viable_points_indices.size() / 2]; + result.points[central_idx].central_enforcer = true; } else { - if (! wrap_around) { - // we can safely use the last enforcer point. - out.push_back(get_center_idx(last_enforcer_start_idx, enforcers_idxs[i])); + size_t central_idx = orig_large_angle_points_indices.size() / 2; + result.points[orig_large_angle_points_indices[central_idx]].central_enforcer = true; + } + } + +} + +// Get index of previous and next perimeter point of the layer. Because SeamCandidates of all polygons of the given layer +// are sequentially stored in the vector, each perimeter contains info about start and end index. These vales are used to +// deduce index of previous and next neigbour in the corresponding perimeter. +std::pair find_previous_and_next_perimeter_point(const std::vector &perimeter_points, + size_t point_index) { + const SeamCandidate ¤t = perimeter_points[point_index]; + int prev = point_index - 1; //for majority of points, it is true that neighbours lie behind and in front of them in the vector + int next = point_index + 1; + + if (point_index == current.perimeter.start_index) { + // if point_index is equal to start, it means that the previous neighbour is at the end + prev = current.perimeter.end_index; + } + + if (point_index == current.perimeter.end_index) { + // if point_index is equal to end, than next neighbour is at the start + next = current.perimeter.start_index; + } + + assert(prev >= 0); + assert(next >= 0); + return {size_t(prev),size_t(next)}; +} + +// Computes all global model info - transforms object, performs raycasting +void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather occlusion meshes: start"; + auto obj_transform = po->trafo_centered(); + indexed_triangle_set triangle_set; + indexed_triangle_set negative_volumes_set; + //add all parts + for (const ModelVolume *model_volume : po->model_object()->volumes) { + if (model_volume->type() == ModelVolumeType::MODEL_PART + || model_volume->type() == ModelVolumeType::NEGATIVE_VOLUME) { + auto model_transformation = model_volume->get_matrix(); + indexed_triangle_set model_its = model_volume->mesh().its; + its_transform(model_its, model_transformation); + if (model_volume->type() == ModelVolumeType::MODEL_PART) { + its_merge(triangle_set, model_its); + } else { + its_merge(negative_volumes_set, model_its); } } } + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather occlusion meshes: end"; - if (wrap_around) { - // Update first center already found. - if (out.empty()) { - // Probably an enforcer around the whole contour. Return nothing. - return out; - } + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: simplify occlusion meshes: start"; - // find last point of the enforcer at the beginning: - size_t idx = 0; - while (enforcers_idxs[idx]+1 == enforcers_idxs[idx+1]) - ++idx; + //simplify raycasting mesh + its_quadric_edge_collapse(triangle_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr, nullptr, + nullptr); + triangle_set = its_subdivide(triangle_set, SeamPlacer::raycasting_subdivision_target_length); - float t_s = lengths[last_enforcer_start_idx]; - float t_e = lengths[idx]; - float half_dist = 0.5f * (t_e + lengths.back() - t_s); - float t_c = (half_dist > t_e) ? t_s + half_dist : t_e - half_dist; + //simplify negative volumes + its_quadric_edge_collapse(negative_volumes_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr, + nullptr, + nullptr); + negative_volumes_set = its_subdivide(negative_volumes_set, SeamPlacer::raycasting_subdivision_target_length); - auto it = std::lower_bound(lengths.begin(), lengths.end(), t_c); - out[0] = it - lengths.begin(); - if (out[0] == lengths.size() - 1) - --out[0]; - assert(out[0] < lengths.size() - 1); - } - return out; + size_t negative_volumes_start_index = triangle_set.indices.size(); + its_merge(triangle_set, negative_volumes_set); + its_transform(triangle_set, obj_transform); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: simplify occlusion meshes: end"; + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer:build AABB tree: start"; + auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices, + triangle_set.indices); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer:build AABB tree: end"; + result.model = triangle_set; + result.model_tree = raycasting_tree; + result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set, negative_volumes_start_index); + +#ifdef DEBUG_FILES + auto filename = debug_out_path(("visiblity_of_" + std::to_string(po->id().id) + ".obj").c_str()); + result.debug_export(triangle_set, filename.c_str()); +#endif } +void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: start"; + auto obj_transform = po->trafo_centered(); -void SeamPlacer::apply_custom_seam(const Polygon& polygon, size_t po_idx, - std::vector& penalties, - const std::vector& lengths, - int layer_id, SeamPosition seam_position) const -{ - if (! is_custom_seam_on_layer(layer_id, po_idx)) - return; + for (const ModelVolume *mv : po->model_object()->volumes) { + if (mv->is_seam_painted()) { + auto model_transformation = obj_transform * mv->get_matrix(); - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs); + indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER); + its_transform(enforcers, model_transformation); + its_merge(result.enforcers, enforcers); - for (size_t i : enforcers_idxs) { - assert(i < penalties.size()); - penalties[i] -= float(ENFORCER_BLOCKER_PENALTY); - } - for (size_t i : blockers_idxs) { - assert(i < penalties.size()); - penalties[i] += float(ENFORCER_BLOCKER_PENALTY); - } - if (seam_position == spAligned) { - std::vector enf_centers = find_enforcer_centers(polygon, lengths, enforcers_idxs); - for (size_t idx : enf_centers) { - assert(idx < penalties.size()); - penalties[idx] += ENFORCER_CENTER_PENALTY; + indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER); + its_transform(blockers, model_transformation); + its_merge(result.blockers, blockers); } } -#if 0 - std::ostringstream os; - os << std::setw(3) << std::setfill('0') << layer_id; - int a = scale_(30.); - SVG svg("custom_seam" + os.str() + ".svg", BoundingBox(Point(-a, -a), Point(a, a))); - if (! m_enforcers[po_idx].empty()) - svg.draw(m_enforcers[po_idx][layer_id].polys, "blue"); - if (! m_blockers[po_idx].empty()) - svg.draw(m_blockers[po_idx][layer_id].polys, "red"); + result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices, + result.enforcers.indices); + result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices, + result.blockers.indices); - if (! blockers_idxs.empty()) { - ; + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end"; +} + +struct SeamComparator { + SeamPosition setup; + + SeamComparator(SeamPosition setup) : + setup(setup) { } + // Standard comparator, must respect the requirements of comparators (e.g. give same result on same inputs) for sorting usage + // should return if a is better seamCandidate than b + bool is_first_better(const SeamCandidate &a, const SeamCandidate &b, const Vec2f &preffered_location = Vec2f { 0.0f, + 0.0f }) const { + if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) { + return a.central_enforcer; + } - size_t min_idx = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); + // Blockers/Enforcers discrimination, top priority + if (a.type != b.type) { + return a.type > b.type; + } - for (size_t i=0; i 0.0f || b.overhang > 0.0f) { + return a.overhang < b.overhang; + } + + // prefer hidden points (more than 1 mm inside) + if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) { + return true; + } + if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) { + return false; + } + + if (setup == SeamPosition::spRear) { + return a.position.y() > b.position.y(); + } + + float distance_penalty_a = 1.0f; + float distance_penalty_b = 1.0f; + if (setup == spNearest) { + distance_penalty_a = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + distance_penalty_b = 1.1f - gauss((b.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + } + + //ranges: [0 - 1] (0 - 1.3] [0.1 - 1.1) + float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) + * compute_angle_penalty(a.local_ccw_angle) + * distance_penalty_a; + float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) + * compute_angle_penalty(b.local_ccw_angle) + * distance_penalty_b; + + return penalty_a < penalty_b; } + + // Comparator used during alignment. If there is close potential aligned point, it is compared to the current + // seam point of the perimeter, to find out if the aligned point is not much worse than the current seam + // Also used by the random seam generator. + bool is_first_not_much_worse(const SeamCandidate &a, const SeamCandidate &b) const { + // Blockers/Enforcers discrimination, top priority + if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) { + // Prefer centers of enforcers. + return a.central_enforcer; + } + + if (a.type == EnforcedBlockedSeamPoint::Enforced) { + return true; + } + + if (a.type == EnforcedBlockedSeamPoint::Blocked) { + return false; + } + + if (a.type != b.type) { + return a.type > b.type; + } + + //avoid overhangs + if (a.overhang > 0.0f || b.overhang > 0.0f) { + return a.overhang < b.overhang; + } + + // prefer hidden points (more than 1 mm inside) + if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) { + return true; + } + if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) { + return false; + } + + if (setup == SeamPosition::spRandom) { + return true; + } + + if (setup == SeamPosition::spRear) { + return a.position.y() > b.position.y(); + } + + //ranges: [0 - 1] (0 - 1.3] ; + float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) + * compute_angle_penalty(a.local_ccw_angle); + float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) + * compute_angle_penalty(b.local_ccw_angle); + + return penalty_a <= penalty_b || penalty_a - penalty_b < SeamPlacer::seam_align_score_tolerance; + } + + bool are_similar(const SeamCandidate &a, const SeamCandidate &b) const { + return is_first_not_much_worse(a, b) && is_first_not_much_worse(b, a); + } + + float compute_angle_penalty(float ccw_angle) const { + // This function is used: + // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x))) + // looks scary, but it is gaussian combined with sigmoid, + // so that concave points have much smaller penalty over convex ones + // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png + return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) + + 1.0f / (2 + std::exp(-ccw_angle)); // sigmoid, which heavily favourizes concave angles + } +}; + +#ifdef DEBUG_FILES +void debug_export_points(const std::vector &layers, + const BoundingBox &bounding_box, std::string object_name, const SeamComparator &comparator) { + for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) { + std::string angles_file_name = debug_out_path( + (object_name + "_angles_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG angles_svg {angles_file_name, bounding_box}; + float min_vis = 0; + float max_vis = min_vis; + + float min_weight = std::numeric_limits::min(); + float max_weight = min_weight; + + for (const SeamCandidate &point : layers[layer_idx].points) { + Vec3i color = value_to_rgbi(-PI, PI, point.local_ccw_angle); + std::string fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + "," + + std::to_string(color.z()) + ")"; + angles_svg.draw(scaled(Vec2f(point.position.head<2>())), fill); + min_vis = std::min(min_vis, point.visibility); + max_vis = std::max(max_vis, point.visibility); + + min_weight = std::min(min_weight, -comparator.compute_angle_penalty(point.local_ccw_angle)); + max_weight = std::max(max_weight, -comparator.compute_angle_penalty(point.local_ccw_angle)); + + } + + std::string visiblity_file_name = debug_out_path( + (object_name + "_visibility_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG visibility_svg {visiblity_file_name, bounding_box}; + std::string weights_file_name = debug_out_path( + (object_name + "_weight_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG weight_svg {weights_file_name, bounding_box}; + std::string overhangs_file_name = debug_out_path( + (object_name + "_overhang_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG overhangs_svg {overhangs_file_name, bounding_box}; + + for (const SeamCandidate &point : layers[layer_idx].points) { + Vec3i color = value_to_rgbi(min_vis, max_vis, point.visibility); + std::string visibility_fill = "rgb(" + std::to_string(color.x()) + "," + std::to_string(color.y()) + "," + + std::to_string(color.z()) + ")"; + visibility_svg.draw(scaled(Vec2f(point.position.head<2>())), visibility_fill); + + Vec3i weight_color = value_to_rgbi(min_weight, max_weight, -comparator.compute_angle_penalty(point.local_ccw_angle)); + std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y()) + + "," + + std::to_string(weight_color.z()) + ")"; + weight_svg.draw(scaled(Vec2f(point.position.head<2>())), weight_fill); + + Vec3i overhang_color = value_to_rgbi(-0.5, 0.5, std::clamp(point.overhang, -0.5f, 0.5f)); + std::string overhang_fill = "rgb(" + std::to_string(overhang_color.x()) + "," + + std::to_string(overhang_color.y()) + + "," + + std::to_string(overhang_color.z()) + ")"; + overhangs_svg.draw(scaled(Vec2f(point.position.head<2>())), overhang_fill); + } + } +} +#endif + +// Pick best seam point based on the given comparator +void pick_seam_point(std::vector &perimeter_points, size_t start_index, + const SeamComparator &comparator) { + size_t end_index = perimeter_points[start_index].perimeter.end_index; + + size_t seam_index = start_index; + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index])) { + seam_index = index; + } + } + perimeter_points[start_index].perimeter.seam_index = seam_index; +} + +size_t pick_nearest_seam_point_index(const std::vector &perimeter_points, size_t start_index, + const Vec2f &preffered_location) { + size_t end_index = perimeter_points[start_index].perimeter.end_index; + SeamComparator comparator { spNearest }; + + size_t seam_index = start_index; + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index], preffered_location)) { + seam_index = index; + } + } + return seam_index; +} + +// picks random seam point uniformly, respecting enforcers blockers and overhang avoidance. +void pick_random_seam_point(const std::vector &perimeter_points, size_t start_index) { + SeamComparator comparator { spRandom }; + + // algorithm keeps a list of viable points and their lengths. If it finds a point + // that is much better than the viable_example_index (e.g. better type, no overhang; see is_first_not_much_worse) + // then it throws away stored lists and starts from start + // in the end, the list should contain points with same type (Enforced > Neutral > Blocked) and also only those which are not + // big overhang. + size_t viable_example_index = start_index; + size_t end_index = perimeter_points[start_index].perimeter.end_index; + struct Viable { + // Candidate seam point index. + size_t index; + float edge_length; + Vec3f edge; + }; + std::vector viables; + + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.are_similar(perimeter_points[index], perimeter_points[viable_example_index])) { + // index ok, push info into viables + Vec3f edge_to_next { perimeter_points[index == end_index ? start_index : index + 1].position + - perimeter_points[index].position }; + float dist_to_next = edge_to_next.norm(); + viables.push_back( { index, dist_to_next, edge_to_next }); + } else if (comparator.is_first_not_much_worse(perimeter_points[viable_example_index], + perimeter_points[index])) { + // index is worse then viable_example_index, skip this point + } else { + // index is better than viable example index, update example, clear gathered info, start again + // clear up all gathered info, start from scratch, update example index + viable_example_index = index; + viables.clear(); + + Vec3f edge_to_next = (perimeter_points[index == end_index ? start_index : index + 1].position + - perimeter_points[index].position); + float dist_to_next = edge_to_next.norm(); + viables.push_back( { index, dist_to_next, edge_to_next }); + } + } + + // now pick random point from the stored options + float len_sum = std::accumulate(viables.begin(), viables.end(), 0.0f, [](const float acc, const Viable &v) { + return acc + v.edge_length; + }); + float picked_len = len_sum * (rand() / (float(RAND_MAX) + 1)); + + size_t point_idx = 0; + while (picked_len - viables[point_idx].edge_length > 0) { + picked_len = picked_len - viables[point_idx].edge_length; + point_idx++; + } + + Perimeter &perimeter = perimeter_points[start_index].perimeter; + perimeter.seam_index = viables[point_idx].index; + perimeter.final_seam_position = perimeter_points[perimeter.seam_index].position + + viables[point_idx].edge.normalized() * picked_len; + perimeter.finalized = true; +} + +struct EdgeGridWrapper { + explicit EdgeGridWrapper(ExPolygons ex_polys) : + ex_polys(ex_polys) { + + grid.create(this->ex_polys, distance_field_resolution); + grid.calculate_sdf(); + } + const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5); + EdgeGrid::Grid grid; + ExPolygons ex_polys; +} +; + +EdgeGridWrapper compute_layer_merged_edge_grid(const Layer *layer) { + static const float eps = float(scale_(layer->object()->config().slice_closing_radius.value)); + // merge with offset + ExPolygons merged = layer->merged(eps); + // ofsset back + ExPolygons layer_outline = offset_ex(merged, -eps); + return EdgeGridWrapper(layer_outline); +} + +} // namespace SeamPlacerImpl + +// Parallel process and extract each perimeter polygon of the given print object. +// Gather SeamCandidates of each layer into vector and build KDtree over them +// Store results in the SeamPlacer variables m_seam_per_object +void SeamPlacer::gather_seam_candidates(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info, const SeamPosition configured_seam_preference) { + using namespace SeamPlacerImpl; + + PrintObjectSeamData &seam_data = m_seam_per_object.emplace(po, PrintObjectSeamData { }).first->second; + seam_data.layers.resize(po->layer_count()); + + tbb::parallel_for(tbb::blocked_range(0, po->layers().size()), + [po, configured_seam_preference, &global_model_info, &seam_data](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + PrintObjectSeamData::LayerSeams &layer_seams = seam_data.layers[layer_idx]; + const Layer *layer = po->get_layer(layer_idx); + auto unscaled_z = layer->slice_z; + std::vector regions; + //NOTE corresponding region ptr may be null, if the layer has zero perimeters + Polygons polygons = extract_perimeter_polygons(layer, configured_seam_preference, regions); + for (size_t poly_index = 0; poly_index < polygons.size(); ++poly_index) { + process_perimeter_polygon(polygons[poly_index], unscaled_z, + regions[poly_index], global_model_info, layer_seams); + } + auto functor = SeamCandidateCoordinateFunctor { layer_seams.points }; + seam_data.layers[layer_idx].points_tree = + std::make_unique(functor, + layer_seams.points.size()); + } + } + ); +} + +void SeamPlacer::calculate_candidates_visibility(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info) { + using namespace SeamPlacerImpl; + + std::vector &layers = m_seam_per_object[po].layers; + tbb::parallel_for(tbb::blocked_range(0, layers.size()), + [&layers, &global_model_info](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + for (auto &perimeter_point : layers[layer_idx].points) { + perimeter_point.visibility = global_model_info.calculate_point_visibility( + perimeter_point.position); + } + } + }); +} + +void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po) { + using namespace SeamPlacerImpl; + + std::vector &layers = m_seam_per_object[po].layers; + tbb::parallel_for(tbb::blocked_range(0, layers.size()), + [po, &layers](tbb::blocked_range r) { + std::unique_ptr prev_layer_grid; + if (r.begin() > 0) { // previous layer exists + prev_layer_grid = std::make_unique( + compute_layer_merged_edge_grid(po->layers()[r.begin() - 1])); + } + + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + bool layer_has_multiple_loops = + layers[layer_idx].points[0].perimeter.end_index + < layers[layer_idx].points.size() - 1; + std::unique_ptr current_layer_grid = std::make_unique( + compute_layer_merged_edge_grid(po->layers()[layer_idx])); + + for (SeamCandidate &perimeter_point : layers[layer_idx].points) { + Point point = Point::new_scale(Vec2f { perimeter_point.position.head<2>() }); + if (prev_layer_grid.get() != nullptr) { + coordf_t overhang_dist; + prev_layer_grid->grid.signed_distance(point, scaled(perimeter_point.perimeter.flow_width), + overhang_dist); + perimeter_point.overhang = + unscale(overhang_dist) - perimeter_point.perimeter.flow_width; + } + + if (layer_has_multiple_loops) { // search for embedded perimeter points (points hidden inside the print ,e.g. multimaterial join, best position for seam) + coordf_t layer_embedded_distance; + current_layer_grid->grid.signed_distance(point, scaled(1.0f), + layer_embedded_distance); + perimeter_point.embedded_distance = unscale(layer_embedded_distance); + } + } + + prev_layer_grid.swap(current_layer_grid); + } + } + ); +} + +// Estimates, if there is good seam point in the layer_idx which is close to last_point_pos +// uses comparator.is_first_not_much_worse method to compare current seam with the closest point +// (if current seam is too far away ) +// If the current chosen stream is close enough, it is stored in seam_string. returns true and updates last_point_pos +// If the closest point is good enough to replace current chosen seam, it is stored in potential_string_seams, returns true and updates last_point_pos +// Otherwise does nothing, returns false +// Used by align_seam_points(). +bool SeamPlacer::find_next_seam_in_layer( + const std::vector &layers, + std::pair &last_point_indexes, + const size_t layer_idx, const float slice_z, + const SeamPlacerImpl::SeamComparator &comparator, + std::vector> &seam_string) const { + using namespace SeamPlacerImpl; + + const SeamCandidate &last_point = layers[last_point_indexes.first].points[last_point_indexes.second]; + + Vec3f projected_position { last_point.position.x(), last_point.position.y(), slice_z }; + std::vector nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position, + SeamPlacer::seam_align_tolerable_dist); + + if (nearby_points_indices.empty()) { + return false; + } + + size_t best_nearby_point_index = nearby_points_indices[0]; + size_t nearest_point_index = nearby_points_indices[0]; + + // Now find best nearby point, nearest point, and corresponding indices + for (const size_t &nearby_point_index : nearby_points_indices) { + const SeamCandidate &point = layers[layer_idx].points[nearby_point_index]; + if (point.perimeter.finalized) { + continue; // skip over finalized perimeters, try to find some that is not finalized + } + if (comparator.is_first_better(point, layers[layer_idx].points[best_nearby_point_index], + projected_position.head<2>()) + || layers[layer_idx].points[best_nearby_point_index].perimeter.finalized) { + best_nearby_point_index = nearby_point_index; + } + if ((point.position - projected_position).squaredNorm() + < (layers[layer_idx].points[nearest_point_index].position - projected_position).squaredNorm() + || layers[layer_idx].points[nearest_point_index].perimeter.finalized) { + nearest_point_index = nearby_point_index; + } + } + + const SeamCandidate &best_nearby_point = layers[layer_idx].points[best_nearby_point_index]; + const SeamCandidate &nearest_point = layers[layer_idx].points[nearest_point_index]; + + if (nearest_point.perimeter.finalized) { + //all points are from already finalized perimeter, skip + return false; + } + + //from the nearest_point, deduce index of seam in the next layer + const SeamCandidate &next_layer_seam = layers[layer_idx].points[nearest_point.perimeter.seam_index]; + + // First try to pick central enforcer if any present + if (next_layer_seam.central_enforcer + && (next_layer_seam.position - projected_position).squaredNorm() + < sqr(3 * SeamPlacer::seam_align_tolerable_dist)) { + last_point_indexes = std::pair { layer_idx, nearest_point.perimeter.seam_index }; + seam_string.push_back(last_point_indexes); + return true; + } + + // Next compare nearest and nearby point. If they are similar pick nearest, Otherwise expect curvy lines on smooth surfaces like chimney of benchy model + // We also compare it to the last point, to detect sharp changes in the scoring - that points to change in the model geometry and string should be ended. + if (comparator.are_similar(nearest_point, best_nearby_point) + && comparator.is_first_not_much_worse(nearest_point, next_layer_seam) + && comparator.are_similar(last_point, nearest_point)) { + last_point_indexes = std::pair { layer_idx, nearest_point_index }; + seam_string.push_back(last_point_indexes); + return true; + } + // If nearest point is not good enough, try it with the best nearby point. + if (comparator.is_first_not_much_worse(best_nearby_point, next_layer_seam) + && comparator.are_similar(last_point, nearest_point)) { + last_point_indexes = std::pair { layer_idx, best_nearby_point_index }; + seam_string.push_back(last_point_indexes); + return true; + } + + return false; +} + +// clusters already chosen seam points into strings across multiple layers, and then +// aligns the strings via polynomial fit +// Does not change the positions of the SeamCandidates themselves, instead stores +// the new aligned position into the shared Perimeter structure of each perimeter +// Note that this position does not necesarilly lay on the perimeter. +void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) { + using namespace SeamPlacerImpl; + + // Prepares Debug files for writing. +#ifdef DEBUG_FILES + Slic3r::CNumericLocalesSetter locales_setter; + auto clusters_f = debug_out_path(("seam_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str()); + FILE *clusters = boost::nowide::fopen(clusters_f.c_str(), "w"); + if (clusters == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << clusters_f << " for writing"; + return; + } + auto aligned_f = debug_out_path(("aligned_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str()); + FILE *aligns = boost::nowide::fopen(aligned_f.c_str(), "w"); + if (aligns == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << clusters_f << " for writing"; + return; + } +#endif + + //gather vector of all seams on the print_object - pair of layer_index and seam__index within that layer + const std::vector &layers = m_seam_per_object[po].layers; + std::vector> seams; + for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) { + const std::vector &layer_perimeter_points = layers[layer_idx].points; + size_t current_point_index = 0; + while (current_point_index < layer_perimeter_points.size()) { + seams.emplace_back(layer_idx, layer_perimeter_points[current_point_index].perimeter.seam_index); + current_point_index = layer_perimeter_points[current_point_index].perimeter.end_index + 1; + } + } + + //sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice + std::sort(seams.begin(), seams.end(), + [&comparator, &layers](const std::pair &left, const std::pair &right) { + return comparator.is_first_better(layers[left.first].points[left.second], + layers[right.first].points[right.second]); + } + ); + + //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment + // Keeping the vectors outside, so with a bit of luck they will not get reallocated after couple of for loop iterations. + std::vector> seam_string; + std::vector observations; + std::vector observation_points; + std::vector weights; + for (const std::pair &seam : seams) { + size_t layer_idx = seam.first; + size_t seam_index = seam.second; + const std::vector &layer_perimeter_points = layers[layer_idx].points; + if (layer_perimeter_points[seam_index].perimeter.finalized) { + // This perimeter is already aligned, skip seam + continue; + } else { + + //initialize searching for seam string - cluster of nearby seams on previous and next layers + int skips = SeamPlacer::seam_align_tolerable_skips / 2; + int next_layer = layer_idx + 1; + std::pair last_point_indexes = std::pair(layer_idx, seam_index); + + seam_string = { std::pair(layer_idx, seam_index) }; + + //find seams or potential seams in forward direction; there is a budget of skips allowed + while (skips >= 0 && next_layer < int(layers.size())) { + if (find_next_seam_in_layer(layers, last_point_indexes, next_layer, + float(po->get_layer(next_layer)->slice_z), comparator, seam_string)) { + //String added, last_point_pos updated, nothing to be done + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer++; + } + + //do additional check in back direction + next_layer = layer_idx - 1; + skips = SeamPlacer::seam_align_tolerable_skips / 2; + last_point_indexes = std::pair(layer_idx, seam_index); + while (skips >= 0 && next_layer >= 0) { + if (find_next_seam_in_layer(layers, last_point_indexes, next_layer, + float(po->get_layer(next_layer)->slice_z), comparator, seam_string)) { + //String added, last_point_pos updated, nothing to be done + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer--; + } + + if (seam_string.size() < seam_align_minimum_string_seams) { + //string NOT long enough to be worth aligning, skip + continue; + } + + // String is long enough, all string seams and potential string seams gathered, now do the alignment + //sort by layer index + std::sort(seam_string.begin(), seam_string.end(), + [](const std::pair &left, const std::pair &right) { + return left.first < right.first; + }); + + // gather all positions of seams and their weights (weights are derived as negative penalty, they are made positive in next step) + observations.resize(seam_string.size()); + observation_points.resize(seam_string.size()); + weights.resize(seam_string.size()); + + //gather points positions and weights + // The algorithm uses only angle to compute penalty, to enforce snapping to sharp corners, if they are present + // after several experiments approach that gives best results is to snap the weight to one for sharp corners, and + // leave it small for others. However, this can result in non-smooth line over area with a lot of unaligned sharp corners. + for (size_t index = 0; index < seam_string.size(); ++index) { + Vec3f pos = layers[seam_string[index].first].points[seam_string[index].second].position; + observations[index] = pos.head<2>(); + observation_points[index] = pos.z(); + weights[index] = + (comparator.compute_angle_penalty( + layers[seam_string[index].first].points[seam_string[index].second].local_ccw_angle) + < comparator.compute_angle_penalty(0.4f * float(PI))) ? 1.0f : 0.1f; + } + + // Curve Fitting + size_t number_of_segments = std::max(size_t(1), + size_t(observations.size() / SeamPlacer::seam_align_seams_per_segment)); + auto curve = Geometry::fit_cubic_bspline(observations, observation_points, weights, number_of_segments); + + // Do alignment - compute fitted point for each point in the string from its Z coord, and store the position into + // Perimeter structure of the point; also set flag aligned to true + for (size_t index = 0; index < seam_string.size(); ++index) { + const auto &pair = seam_string[index]; + const float t = weights[index]; + Vec3f current_pos = layers[pair.first].points[pair.second].position; + Vec2f fitted_pos = curve.get_fitted_value(current_pos.z()); + + //interpolate between current and fitted position, prefer current pos for large weights. + Vec3f final_position = t * current_pos + (1 - t) * to_3d(fitted_pos, current_pos.z()); + + Perimeter &perimeter = layers[pair.first].points[pair.second].perimeter; + perimeter.seam_index = pair.second; + perimeter.final_seam_position = final_position; + perimeter.finalized = true; + } + +#ifdef DEBUG_FILES + auto randf = []() { + return float(rand()) / float(RAND_MAX); + }; + Vec3f color { randf(), randf(), randf() }; + for (size_t i = 0; i < seam_string.size(); ++i) { + auto orig_seam = layers[seam_string[i].first].points[seam_string[i].second]; + fprintf(clusters, "v %f %f %f %f %f %f \n", orig_seam.position[0], + orig_seam.position[1], + orig_seam.position[2], color[0], color[1], + color[2]); + } + + color = Vec3f { randf(), randf(), randf() }; + for (size_t i = 0; i < seam_string.size(); ++i) { + const Perimeter &perimeter = layers[seam_string[i].first].points[seam_string[i].second].perimeter; + fprintf(aligns, "v %f %f %f %f %f %f \n", perimeter.final_seam_position[0], + perimeter.final_seam_position[1], + perimeter.final_seam_position[2], color[0], color[1], + color[2]); + } +#endif + } + } + +#ifdef DEBUG_FILES + fclose(clusters); + fclose(aligns); #endif } +void SeamPlacer::init(const Print &print) { + using namespace SeamPlacerImpl; + m_seam_per_object.clear(); + for (const PrintObject *po : print.objects()) { -std::optional SeamHistory::get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb) -{ - assert(layer_id >= m_layer_id || layer_id == 0); - if (layer_id != m_layer_id) { - // Get seam was called for different layer than last time. - if (layer_id == 0) // seq printing - m_data_this_layer.clear(); - m_data_last_layer = m_data_this_layer; - m_data_this_layer.clear(); - m_layer_id = layer_id; - } + SeamPosition configured_seam_preference = po->config().seam_position.value; + SeamComparator comparator { configured_seam_preference }; - std::optional out; + GlobalModelInfo global_model_info { }; + gather_enforcers_blockers(global_model_info, po); - auto seams_it = m_data_last_layer.find(po); - if (seams_it == m_data_last_layer.end()) - return out; - - const std::vector& seam_data_po = seams_it->second; - - // Find a bounding-box on the last layer that is close to one we see now. - double min_score = std::numeric_limits::max(); - for (const SeamPoint& sp : seam_data_po) { - const BoundingBox& bb = sp.m_island_bb; - - if (! bb.overlap(island_bb)) { - // This bb does not even overlap. It is likely unrelated. - continue; + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { + compute_global_occlusion(global_model_info, po); } - double score = std::pow(bb.min(0) - island_bb.min(0), 2.) - + std::pow(bb.min(1) - island_bb.min(1), 2.) - + std::pow(bb.max(0) - island_bb.max(0), 2.) - + std::pow(bb.max(1) - island_bb.max(1), 2.); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather_seam_candidates: start"; + gather_seam_candidates(po, global_model_info, configured_seam_preference); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather_seam_candidates: end"; - if (score < min_score) { - min_score = score; - out = sp.m_pos; + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_candidates_visibility : start"; + calculate_candidates_visibility(po, global_model_info); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_candidates_visibility : end"; } + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_overhangs and layer embdedding : start"; + calculate_overhangs_and_layer_embedding(po); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_overhangs and layer embdedding: end"; + + if (configured_seam_preference != spNearest) { // For spNearest, the seam is picked in the place_seam method with actual nozzle position information + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: pick_seam_point : start"; + //pick seam point + std::vector &layers = m_seam_per_object[po].layers; + tbb::parallel_for(tbb::blocked_range(0, layers.size()), + [&layers, configured_seam_preference, comparator](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + std::vector &layer_perimeter_points = layers[layer_idx].points; + for (size_t current = 0; current < layer_perimeter_points.size(); + current = layer_perimeter_points[current].perimeter.end_index + 1) + if (configured_seam_preference == spRandom) + pick_random_seam_point(layer_perimeter_points, current); + else + pick_seam_point(layer_perimeter_points, current, comparator); + } + }); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: pick_seam_point : end"; + } + + if (configured_seam_preference == spAligned) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: align_seam_points : start"; + align_seam_points(po, comparator); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: align_seam_points : end"; + } + +#ifdef DEBUG_FILES + debug_export_points(layers, po->bounding_box(), std::to_string(po->id().id), + comparator); +#endif + } +} + +void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, + const Point &last_pos) const { + using namespace SeamPlacerImpl; + const PrintObject *po = layer->object(); + // Must not be called with supprot layer. + assert(dynamic_cast(layer) == nullptr); + // Object layer IDs are incremented by the number of raft layers. + assert(layer->id() >= po->slicing_parameters().raft_layers()); + const size_t layer_index = layer->id() - po->slicing_parameters().raft_layers(); + const double unscaled_z = layer->slice_z; + + const PrintObjectSeamData::LayerSeams &layer_perimeters = + m_seam_per_object.find(layer->object())->second.layers[layer_index]; + + // Find the closest perimeter in the SeamPlacer to the first point of this loop. + size_t closest_perimeter_point_index; + { + const Point &fp = loop.first_point(); + Vec2f unscaled_p = unscaled(fp); + closest_perimeter_point_index = find_closest_point(*layer_perimeters.points_tree.get(), + to_3d(unscaled_p, float(unscaled_z))); } - return out; -} + Vec3f seam_position; + size_t seam_index; + if (const Perimeter &perimeter = layer_perimeters.points[closest_perimeter_point_index].perimeter; + perimeter.finalized) { + seam_position = perimeter.final_seam_position; + seam_index = perimeter.seam_index; + } else { + seam_index = + po->config().seam_position == spNearest ? + pick_nearest_seam_point_index(layer_perimeters.points, perimeter.start_index, + unscaled(last_pos)) : + perimeter.seam_index; + seam_position = layer_perimeters.points[seam_index].position; + } + Point seam_point = Point::new_scale(seam_position.x(), seam_position.y()); + if (const SeamCandidate &perimeter_point = layer_perimeters.points[seam_index]; + (po->config().seam_position == spNearest || po->config().seam_position == spAligned) && + loop.role() == ExtrusionRole::erPerimeter && //Hopefully internal perimeter + (seam_position - perimeter_point.position).squaredNorm() < 4.0f && // seam is on perimeter point + perimeter_point.local_ccw_angle < -EPSILON // In concave angles + ) { // In this case, we are at internal perimeter, where the external perimeter has seam in concave angle. We want to align + // the internal seam into the concave corner, and not on the perpendicular projection on the closest edge (which is what the split_at function does) + size_t index_of_prev = + seam_index == perimeter_point.perimeter.start_index ? + perimeter_point.perimeter.end_index : + seam_index - 1; + size_t index_of_next = + seam_index == perimeter_point.perimeter.end_index ? + perimeter_point.perimeter.start_index : + seam_index + 1; -void SeamHistory::add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb) -{ - m_data_this_layer[po].push_back({pos, island_bb});; -} + Vec2f dir_to_middle = + ((perimeter_point.position - layer_perimeters.points[index_of_prev].position).head<2>().normalized() + + (perimeter_point.position - layer_perimeters.points[index_of_next].position).head<2>().normalized()) + * 0.5; + auto [_, projected_point] = loop.get_closest_path_and_point(seam_point, true); + //get closest projected point, determine depth of the seam point. + float depth = (float) unscale(Point(seam_point - projected_point)).norm(); + float angle_factor = cos(-perimeter_point.local_ccw_angle / 2.0f); // There are some nice geometric identities in determination of the correct depth of new seam point. + //overshoot the target depth, in concave angles it will correctly snap to the corner; TODO: find out why such big overshoot is needed. + Vec2f final_pos = perimeter_point.position.head<2>() + (1.4142 * depth / angle_factor) * dir_to_middle; + seam_point = Point::new_scale(final_pos.x(), final_pos.y()); + } - -void SeamHistory::clear() -{ - m_layer_id = 0; - m_data_last_layer.clear(); - m_data_this_layer.clear(); -} - + if (!loop.split_at_vertex(seam_point)) { + // The point is not in the original loop. + // Insert it. + loop.split_at(seam_point, true); + } } + +} // namespace Slic3r diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 0d72939b2..70881d558 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -3,12 +3,16 @@ #include #include +#include +#include +#include "libslic3r/libslic3r.h" #include "libslic3r/ExtrusionEntity.hpp" #include "libslic3r/Polygon.hpp" #include "libslic3r/PrintConfig.hpp" #include "libslic3r/BoundingBox.hpp" #include "libslic3r/AABBTreeIndirect.hpp" +#include "libslic3r/KDTreeIndirect.hpp" namespace Slic3r { @@ -16,119 +20,148 @@ class PrintObject; class ExtrusionLoop; class Print; class Layer; -namespace EdgeGrid { class Grid; } +namespace EdgeGrid { +class Grid; +} -class SeamHistory { -public: - SeamHistory() { clear(); } - std::optional get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb); - void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb); - void clear(); +namespace SeamPlacerImpl { -private: - struct SeamPoint { - Point m_pos; - BoundingBox m_island_bb; - }; +struct GlobalModelInfo; +struct SeamComparator; - std::map> m_data_last_layer; - std::map> m_data_this_layer; - size_t m_layer_id; +enum class EnforcedBlockedSeamPoint { + Blocked = 0, + Neutral = 1, + Enforced = 2, }; +// struct representing single perimeter loop +struct Perimeter { + size_t start_index; + size_t end_index; //inclusive! + size_t seam_index; + float flow_width; + // During alignment, a final position may be stored here. In that case, finalized is set to true. + // Note that final seam position is not limited to points of the perimeter loop. In theory it can be any position + // Random position also uses this flexibility to set final seam point position + bool finalized = false; + Vec3f final_seam_position; +}; + +//Struct over which all processing of perimeters is done. For each perimeter point, its respective candidate is created, +// then all the needed attributes are computed and finally, for each perimeter one point is chosen as seam. +// This seam position can be then further aligned +struct SeamCandidate { + SeamCandidate(const Vec3f &pos, Perimeter &perimeter, + float local_ccw_angle, + EnforcedBlockedSeamPoint type) : + position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), embedded_distance(0.0f), local_ccw_angle( + local_ccw_angle), type(type), central_enforcer(false) { + } + const Vec3f position; + // pointer to Perimeter loop of this point. It is shared across all points of the loop + Perimeter &perimeter; + float visibility; + float overhang; + // distance inside the merged layer regions, for detecting perimeter points which are hidden indside the print (e.g. multimaterial join) + // Negative sign means inside the print, comes from EdgeGrid structure + float embedded_distance; + float local_ccw_angle; + EnforcedBlockedSeamPoint type; + bool central_enforcer; //marks this candidate as central point of enforced segment on the perimeter - important for alignment +}; + +struct FaceVisibilityInfo { + float visibility; +}; + +struct SeamCandidateCoordinateFunctor { + SeamCandidateCoordinateFunctor(const std::vector &seam_candidates) : + seam_candidates(seam_candidates) { + } + const std::vector &seam_candidates; + float operator()(size_t index, size_t dim) const { + return seam_candidates[index].position[dim]; + } +}; +} // namespace SeamPlacerImpl + +struct PrintObjectSeamData +{ + using SeamCandidatesTree = KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>; + + struct LayerSeams + { + Slic3r::deque perimeters; + std::vector points; + std::unique_ptr points_tree; + }; + // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter + std::vector layers; + // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD + // tree of all points of the given layer + + void clear() + { + layers.clear(); + } +}; class SeamPlacer { public: - void init(const Print& print); + static constexpr size_t raycasting_decimation_target_triangle_count = 10000; + static constexpr float raycasting_subdivision_target_length = 2.0f; + //square of number of rays per triangle + static constexpr size_t sqr_rays_per_triangle = 7; - // When perimeters are printed, first call this function with the respective - // external perimeter. SeamPlacer will find a location for its seam and remember it. - // Subsequent calls to get_seam will return this position. + // arm length used during angles computation + static constexpr float polygon_local_angles_arm_distance = 0.5f; + // increases angle importance at the cost of deacreasing visibility info importance. must be > 0 + static constexpr float additional_angle_importance = 0.6f; - void plan_perimeters(const std::vector perimeters, - const Layer& layer, SeamPosition seam_position, - Point last_pos, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid); + // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer + static constexpr float enforcer_blocker_distance_tolerance = 0.35f; + // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size + static constexpr float enforcer_oversampling_distance = 0.2f; - void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter, - const EdgeGrid::Grid* lower_layer_edge_grid); - + // When searching for seam clusters for alignment: + // following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer + static constexpr float seam_align_score_tolerance = 0.5f; + // seam_align_tolerable_dist - if next layer closes point is too far away, break string + static constexpr float seam_align_tolerable_dist = 1.0f; + // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. + // this param limits the number of allowed skips + static constexpr size_t seam_align_tolerable_skips = 4; + // minimum number of seams needed in cluster to make alignment happen + static constexpr size_t seam_align_minimum_string_seams = 6; + // points covered by spline; determines number of splines for the given string + static constexpr size_t seam_align_seams_per_segment = 8; - using TreeType = AABBTreeIndirect::Tree<2, coord_t>; - using AlignedBoxType = Eigen::AlignedBox; + //The following data structures hold all perimeter points for all PrintObject. + std::unordered_map m_seam_per_object; + + void init(const Print &print); + + void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const; private: - - // When given an external perimeter (!), returns the seam. - Point calculate_seam(const Layer& layer, const SeamPosition seam_position, - const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest); - - struct CustomTrianglesPerLayer { - Polygons polys; - TreeType tree; - }; - - // Just a cache to save some lookups. - const Layer* m_last_layer_po = nullptr; - coordf_t m_last_print_z = -1.; - const PrintObject* m_last_po = nullptr; - - struct SeamPoint { - Point pt; - bool precalculated = false; - bool external = false; - const Layer* layer = nullptr; - SeamPosition seam_position; - const PrintObject* po = nullptr; - }; - std::vector m_plan; - size_t m_plan_idx; - - std::vector> m_enforcers; - std::vector> m_blockers; - std::vector m_po_list; - - //std::map m_last_seam_position; - SeamHistory m_seam_history; - - // Get indices of points inside enforcers and blockers. - void get_enforcers_and_blockers(size_t layer_id, - const Polygon& polygon, - size_t po_id, - std::vector& enforcers_idxs, - std::vector& blockers_idxs) const; - - // Apply penalties to points inside enforcers/blockers. - void apply_custom_seam(const Polygon& polygon, size_t po_id, - std::vector& penalties, - const std::vector& lengths, - int layer_id, SeamPosition seam_position) const; - - // Return random point of a polygon. The distribution will be uniform - // along the contour and account for enforcers and blockers. - Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id, - bool* saw_custom = nullptr) const; - - // Is there any enforcer/blocker on this layer? - bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const { - return is_custom_enforcer_on_layer(layer_id, po_idx) - || is_custom_blocker_on_layer(layer_id, po_idx); - } - - bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const { - return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty()); - } - - bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const { - return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty()); - } + void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info, + const SeamPosition configured_seam_preference); + void calculate_candidates_visibility(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info); + void calculate_overhangs_and_layer_embedding(const PrintObject *po); + void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator); + bool find_next_seam_in_layer( + const std::vector &layers, + std::pair &last_point_indexes, + const size_t layer_idx, const float slice_z, + const SeamPlacerImpl::SeamComparator &comparator, + std::vector> &seam_string) const; }; - -} +} // namespace Slic3r #endif // libslic3r_SeamPlacer_hpp_ diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp index c825f8254..82ffbd8d1 100644 --- a/src/libslic3r/Geometry.hpp +++ b/src/libslic3r/Geometry.hpp @@ -36,9 +36,9 @@ enum Orientation static inline Orientation orient(const Point &a, const Point &b, const Point &c) { static_assert(sizeof(coord_t) * 2 == sizeof(int64_t), "orient works with 32 bit coordinates"); - int64_t u = int64_t(b(0)) * int64_t(c(1)) - int64_t(b(1)) * int64_t(c(0)); - int64_t v = int64_t(a(0)) * int64_t(c(1)) - int64_t(a(1)) * int64_t(c(0)); - int64_t w = int64_t(a(0)) * int64_t(b(1)) - int64_t(a(1)) * int64_t(b(0)); + int64_t u = int64_t(b.x()) * int64_t(c.y()) - int64_t(b.y()) * int64_t(c.x()); + int64_t v = int64_t(a.x()) * int64_t(c.y()) - int64_t(a.y()) * int64_t(c.x()); + int64_t w = int64_t(a.x()) * int64_t(b.y()) - int64_t(a.y()) * int64_t(b.x()); int64_t d = u - v + w; return (d > 0) ? ORIENTATION_CCW : ((d == 0) ? ORIENTATION_COLINEAR : ORIENTATION_CW); } diff --git a/src/libslic3r/Geometry/Bicubic.hpp b/src/libslic3r/Geometry/Bicubic.hpp new file mode 100644 index 000000000..98c6f8bb2 --- /dev/null +++ b/src/libslic3r/Geometry/Bicubic.hpp @@ -0,0 +1,291 @@ +#ifndef BICUBIC_HPP +#define BICUBIC_HPP + +#include +#include +#include + +#include + +namespace Slic3r { + +namespace Geometry { + +namespace BicubicInternal { +// Linear kernel, to be able to test cubic methods with hat kernels. +template +struct LinearKernel +{ + typedef T FloatType; + + static T a00() { + return T(0.); + } + static T a01() { + return T(0.); + } + static T a02() { + return T(0.); + } + static T a03() { + return T(0.); + } + static T a10() { + return T(1.); + } + static T a11() { + return T(-1.); + } + static T a12() { + return T(0.); + } + static T a13() { + return T(0.); + } + static T a20() { + return T(0.); + } + static T a21() { + return T(1.); + } + static T a22() { + return T(0.); + } + static T a23() { + return T(0.); + } + static T a30() { + return T(0.); + } + static T a31() { + return T(0.); + } + static T a32() { + return T(0.); + } + static T a33() { + return T(0.); + } +}; + +// Interpolation kernel aka Catmul-Rom aka Keyes kernel. +template +struct CubicCatmulRomKernel +{ + typedef T FloatType; + + static T a00() { + return 0; + } + static T a01() { + return T( -0.5); + } + static T a02() { + return T( 1.); + } + static T a03() { + return T( -0.5); + } + static T a10() { + return T( 1.); + } + static T a11() { + return 0; + } + static T a12() { + return T( -5. / 2.); + } + static T a13() { + return T( 3. / 2.); + } + static T a20() { + return 0; + } + static T a21() { + return T( 0.5); + } + static T a22() { + return T( 2.); + } + static T a23() { + return T( -3. / 2.); + } + static T a30() { + return 0; + } + static T a31() { + return 0; + } + static T a32() { + return T( -0.5); + } + static T a33() { + return T( 0.5); + } +}; + +// B-spline kernel +template +struct CubicBSplineKernel +{ + typedef T FloatType; + + static T a00() { + return T( 1. / 6.); + } + static T a01() { + return T( -3. / 6.); + } + static T a02() { + return T( 3. / 6.); + } + static T a03() { + return T( -1. / 6.); + } + static T a10() { + return T( 4. / 6.); + } + static T a11() { + return 0; + } + static T a12() { + return T( -6. / 6.); + } + static T a13() { + return T( 3. / 6.); + } + static T a20() { + return T( 1. / 6.); + } + static T a21() { + return T( 3. / 6.); + } + static T a22() { + return T( 3. / 6.); + } + static T a23() { + return T( -3. / 6.); + } + static T a30() { + return 0; + } + static T a31() { + return 0; + } + static T a32() { + return 0; + } + static T a33() { + return T( 1. / 6.); + } +}; + +template +inline T clamp(T a, T lower, T upper) + { + return (a < lower) ? lower : + (a > upper) ? upper : a; +} +} + +template +struct CubicKernelWrapper +{ + typedef typename Kernel::FloatType FloatType; + + static constexpr size_t kernel_span = 4; + + static FloatType kernel(FloatType x) + { + x = fabs(x); + if (x >= (FloatType) 2.) + return 0.0f; + if (x <= (FloatType) 1.) { + FloatType x2 = x * x; + FloatType x3 = x2 * x; + return Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3; + } + assert(x > (FloatType )1. && x < (FloatType )2.); + x -= (FloatType) 1.; + FloatType x2 = x * x; + FloatType x3 = x2 * x; + return Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3; + } + + static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x) + { + const FloatType x2 = x * x; + const FloatType x3 = x * x * x; + return f0 * (Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3) + + f1 * (Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3) + + f2 * (Kernel::a20() + Kernel::a21() * x + Kernel::a22() * x2 + Kernel::a23() * x3) + + f3 * (Kernel::a30() + Kernel::a31() * x + Kernel::a32() * x2 + Kernel::a33() * x3); + } +}; + +// Linear splines +template +using LinearKernel = CubicKernelWrapper>; + +// Catmul-Rom splines +template +using CubicCatmulRomKernel = CubicKernelWrapper>; + +// Cubic B-splines +template +using CubicBSplineKernel = CubicKernelWrapper>; + +template +static typename KernelWrapper::FloatType cubic_interpolate(const Eigen::ArrayBase &F, + const typename KernelWrapper::FloatType pt) { + typedef typename KernelWrapper::FloatType T; + const int w = int(F.size()); + const int ix = (int) floor(pt); + const T s = pt - T( ix); + + if (ix > 1 && ix + 2 < w) { + // Inside the fully interpolated region. + return KernelWrapper::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s); + } + // Transition region. Extend with a constant function. + auto f = [&F, w](T x) { + return F[BicubicInternal::clamp(x, 0, w - 1)]; + }; + return KernelWrapper::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s); +} + +template +static float bicubic_interpolate(const Eigen::MatrixBase &F, + const Eigen::Matrix &pt) { + typedef typename Kernel::FloatType T; + const int w = F.cols(); + const int h = F.rows(); + const int ix = (int) floor(pt[0]); + const int iy = (int) floor(pt[1]); + const T s = pt[0] - T( ix); + const T t = pt[1] - T( iy); + + if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) { + // Inside the fully interpolated region. + return Kernel::interpolate( + Kernel::interpolate(F(ix - 1, iy - 1), F(ix, iy - 1), F(ix + 1, iy - 1), F(ix + 2, iy - 1), s), + Kernel::interpolate(F(ix - 1, iy), F(ix, iy), F(ix + 1, iy), F(ix + 2, iy), s), + Kernel::interpolate(F(ix - 1, iy + 1), F(ix, iy + 1), F(ix + 1, iy + 1), F(ix + 2, iy + 1), s), + Kernel::interpolate(F(ix - 1, iy + 2), F(ix, iy + 2), F(ix + 1, iy + 2), F(ix + 2, iy + 2), s), t); + } + // Transition region. Extend with a constant function. + auto f = [&F, w, h](int x, int y) { + return F(BicubicInternal::clamp(x, 0, w - 1), BicubicInternal::clamp(y, 0, h - 1)); + }; + return Kernel::interpolate( + Kernel::interpolate(f(ix - 1, iy - 1), f(ix, iy - 1), f(ix + 1, iy - 1), f(ix + 2, iy - 1), s), + Kernel::interpolate(f(ix - 1, iy), f(ix, iy), f(ix + 1, iy), f(ix + 2, iy), s), + Kernel::interpolate(f(ix - 1, iy + 1), f(ix, iy + 1), f(ix + 1, iy + 1), f(ix + 2, iy + 1), s), + Kernel::interpolate(f(ix - 1, iy + 2), f(ix, iy + 2), f(ix + 1, iy + 2), f(ix + 2, iy + 2), s), t); +} + +} //namespace Geometry + +} // namespace Slic3r + +#endif /* BICUBIC_HPP */ diff --git a/src/libslic3r/Geometry/ConvexHull.cpp b/src/libslic3r/Geometry/ConvexHull.cpp index b1ff77f80..2e92535f2 100644 --- a/src/libslic3r/Geometry/ConvexHull.cpp +++ b/src/libslic3r/Geometry/ConvexHull.cpp @@ -1,6 +1,7 @@ #include "libslic3r.h" #include "ConvexHull.hpp" #include "BoundingBox.hpp" +#include "../Geometry.hpp" #include @@ -19,13 +20,13 @@ Polygon convex_hull(Points pts) hull.points.resize(2 * n); // Build lower hull for (int i = 0; i < n; ++ i) { - while (k >= 2 && pts[i].ccw(hull[k-2], hull[k-1]) <= 0) + while (k >= 2 && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW) -- k; hull[k ++] = pts[i]; } // Build upper hull for (int i = n-2, t = k+1; i >= 0; i--) { - while (k >= t && pts[i].ccw(hull[k-2], hull[k-1]) <= 0) + while (k >= t && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW) -- k; hull[k ++] = pts[i]; } @@ -58,7 +59,7 @@ Pointf3s convex_hull(Pointf3s points) Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1)); Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1)); - if (p.ccw(k2, k1) <= 0) + if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW) --k; else break; @@ -76,7 +77,7 @@ Pointf3s convex_hull(Pointf3s points) Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1)); Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1)); - if (p.ccw(k2, k1) <= 0) + if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW) --k; else break; diff --git a/src/libslic3r/Geometry/Curves.hpp b/src/libslic3r/Geometry/Curves.hpp new file mode 100644 index 000000000..6542cb706 --- /dev/null +++ b/src/libslic3r/Geometry/Curves.hpp @@ -0,0 +1,205 @@ +#ifndef SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ +#define SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ + +#include "libslic3r/Point.hpp" +#include "Bicubic.hpp" + +#include + +//#define LSQR_DEBUG + +namespace Slic3r { +namespace Geometry { + +template +struct PolynomialCurve { + Eigen::MatrixXf coefficients; + + Vec3f get_fitted_value(const NumberType value) const { + auto result = Vec::Zero(); + size_t order = this->coefficients.rows() - 1; + auto x = NumberType(1.); + for (size_t index = 0; index < order + 1; ++index, x *= value) + result += x * this->coefficients.col(index); + return result; + } +}; + +//https://towardsdatascience.com/least-square-polynomial-CURVES-using-c-eigen-package-c0673728bd01 +template +PolynomialCurve fit_polynomial(const std::vector> &observations, + const std::vector &observation_points, + const std::vector &weights, size_t order) { + // check to make sure inputs are correct + size_t cols = order + 1; + assert(observation_points.size() >= cols); + assert(observation_points.size() == weights.size()); + assert(observations.size() == weights.size()); + + Eigen::MatrixXf data_points(Dimension, observations.size()); + Eigen::MatrixXf T(observations.size(), cols); + for (size_t i = 0; i < weights.size(); ++i) { + auto squared_weight = sqrt(weights[i]); + data_points.col(i) = observations[i] * squared_weight; + // Populate the matrix + auto x = squared_weight; + auto c = observation_points[i]; + for (size_t j = 0; j < cols; ++j, x *= c) + T(i, j) = x; + } + + const auto QR = T.householderQr(); + Eigen::MatrixXf coefficients(Dimension, cols); + // Solve for linear least square fit + for (size_t dim = 0; dim < Dimension; ++dim) { + coefficients.row(dim) = QR.solve(data_points.row(dim).transpose()); + } + + return {std::move(coefficients)}; +} + +template +struct PiecewiseFittedCurve { + using Kernel = KernelType; + + Eigen::MatrixXf coefficients; + NumberType start; + NumberType segment_size; + size_t endpoints_level_of_freedom; + + Vec get_fitted_value(const NumberType &observation_point) const { + Vec result = Vec::Zero(); + + //find corresponding segment index; expects kernels to be centered + int middle_right_segment_index = floor((observation_point - start) / segment_size); + //find index of first segment that is affected by the point i; this can be deduced from kernel_span + int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1; + for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span); + segment_index++) { + NumberType segment_start = start + segment_index * segment_size; + NumberType normalized_segment_distance = (segment_start - observation_point) / segment_size; + + int parameter_index = segment_index + endpoints_level_of_freedom; + parameter_index = std::clamp(parameter_index, 0, int(coefficients.cols()) - 1); + result += Kernel::kernel(normalized_segment_distance) * coefficients.col(parameter_index); + } + return result; + } +}; + +// observations: data to be fitted by the curve +// observation points: growing sequence of points where the observations were made. +// In other words, for function f(x) = y, observations are y0...yn, and observation points are x0...xn +// weights: how important the observation is +// segments_count: number of segments inside the valid length of the curve +// endpoints_level_of_freedom: number of additional parameters at each end; reasonable values depend on the kernel span +template +PiecewiseFittedCurve fit_curve( + const std::vector> &observations, + const std::vector &observation_points, + const std::vector &weights, + size_t segments_count, + size_t endpoints_level_of_freedom) { + + // check to make sure inputs are correct + assert(segments_count > 0); + assert(observations.size() > 0); + assert(observation_points.size() == observations.size()); + assert(observation_points.size() == weights.size()); + assert(segments_count <= observations.size()); + + //prepare sqrt of weights, which will then be applied to both matrix T and observed data: https://en.wikipedia.org/wiki/Weighted_least_squares + std::vector sqrt_weights(weights.size()); + for (size_t index = 0; index < weights.size(); ++index) { + assert(weights[index] > 0); + sqrt_weights[index] = sqrt(weights[index]); + } + + // prepare result and compute metadata + PiecewiseFittedCurve result { }; + + NumberType valid_length = observation_points.back() - observation_points.front(); + NumberType segment_size = valid_length / NumberType(segments_count); + result.start = observation_points.front(); + result.segment_size = segment_size; + result.endpoints_level_of_freedom = endpoints_level_of_freedom; + + // prepare observed data + // Eigen defaults to column major memory layout. + Eigen::MatrixXf data_points(Dimension, observations.size()); + for (size_t index = 0; index < observations.size(); ++index) { + data_points.col(index) = observations[index] * sqrt_weights[index]; + } + // parameters count is always increased by one to make the parametric space of the curve symmetric. + // without this fix, the end of the curve is less flexible than the beginning + size_t parameters_count = segments_count + 1 + 2 * endpoints_level_of_freedom; + //Create weight matrix T for each point and each segment; + Eigen::MatrixXf T(observation_points.size(), parameters_count); + T.setZero(); + //Fill the weight matrix + for (size_t i = 0; i < observation_points.size(); ++i) { + NumberType observation_point = observation_points[i]; + //find corresponding segment index; expects kernels to be centered + int middle_right_segment_index = floor((observation_point - result.start) / result.segment_size); + //find index of first segment that is affected by the point i; this can be deduced from kernel_span + int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1; + for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span); + segment_index++) { + NumberType segment_start = result.start + segment_index * result.segment_size; + NumberType normalized_segment_distance = (segment_start - observation_point) / result.segment_size; + + int parameter_index = segment_index + endpoints_level_of_freedom; + parameter_index = std::clamp(parameter_index, 0, int(parameters_count) - 1); + T(i, parameter_index) += Kernel::kernel(normalized_segment_distance) * sqrt_weights[i]; + } + } + +#ifdef LSQR_DEBUG + std::cout << "weight matrix: " << std::endl; + for (int obs = 0; obs < observation_points.size(); ++obs) { + std::cout << std::endl; + for (int segment = 0; segment < parameters_count; ++segment) { + std::cout << T(obs, segment) << " "; + } + } + std::cout << std::endl; +#endif + + // Solve for linear least square fit + result.coefficients.resize(Dimension, parameters_count); + const auto QR = T.fullPivHouseholderQr(); + for (size_t dim = 0; dim < Dimension; ++dim) { + result.coefficients.row(dim) = QR.solve(data_points.row(dim).transpose()); + } + + return result; +} + +template +PiecewiseFittedCurve> +fit_cubic_bspline( + const std::vector> &observations, + std::vector observation_points, + std::vector weights, + size_t segments_count, + size_t endpoints_level_of_freedom = 0) { + return fit_curve>(observations, observation_points, weights, segments_count, + endpoints_level_of_freedom); +} + +template +PiecewiseFittedCurve> +fit_catmul_rom_spline( + const std::vector> &observations, + std::vector observation_points, + std::vector weights, + size_t segments_count, + size_t endpoints_level_of_freedom = 0) { + return fit_curve>(observations, observation_points, weights, segments_count, + endpoints_level_of_freedom); +} + +} +} + +#endif /* SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ */ diff --git a/src/libslic3r/KDTreeIndirect.hpp b/src/libslic3r/KDTreeIndirect.hpp index 12e462569..36a157456 100644 --- a/src/libslic3r/KDTreeIndirect.hpp +++ b/src/libslic3r/KDTreeIndirect.hpp @@ -59,7 +59,7 @@ public: CONTINUE_RIGHT = 2, STOP = 4, }; - template + template unsigned int descent_mask(const CoordType &point_coord, const CoordType &search_radius, size_t idx, size_t dimension) const { CoordType dist = point_coord - this->coordinate(idx, dimension); @@ -110,8 +110,8 @@ private: // Partition the input m_nodes at "k" and "dimension" using the QuickSelect method: // https://en.wikipedia.org/wiki/Quickselect - // Items left of the k'th item are lower than the k'th item in the "dimension", - // items right of the k'th item are higher than the k'th item in the "dimension", + // Items left of the k'th item are lower than the k'th item in the "dimension", + // items right of the k'th item are higher than the k'th item in the "dimension", void partition_input(std::vector &input, const size_t dimension, size_t left, size_t right, const size_t k) const { while (left < right) { @@ -231,6 +231,53 @@ size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& poi return find_closest_point(kdtree, point, [](size_t) { return true; }); } +// Find nearby points (spherical neighbourhood) using Euclidian metrics. +template +std::vector find_nearby_points(const KDTreeIndirectType &kdtree, const PointType ¢er, + const typename KDTreeIndirectType::CoordType& max_distance, FilterFn filter) + { + using CoordType = typename KDTreeIndirectType::CoordType; + + struct Visitor { + const KDTreeIndirectType &kdtree; + const PointType center; + const CoordType max_distance_squared; + const FilterFn filter; + std::vector result; + + Visitor(const KDTreeIndirectType &kdtree, const PointType& center, const CoordType &max_distance, + FilterFn filter) : + kdtree(kdtree), center(center), max_distance_squared(max_distance*max_distance), filter(filter) { + } + unsigned int operator()(size_t idx, size_t dimension) { + if (this->filter(idx)) { + auto dist = CoordType(0); + for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) { + CoordType d = center[i] - kdtree.coordinate(idx, i); + dist += d * d; + } + if (dist < max_distance_squared) { + result.push_back(idx); + } + } + return kdtree.descent_mask(center[dimension], max_distance_squared, idx, dimension); + } + } visitor(kdtree, center, max_distance, filter); + + kdtree.visit(visitor); + return visitor.result; +} + +template +std::vector find_nearby_points(const KDTreeIndirectType &kdtree, const PointType ¢er, + const typename KDTreeIndirectType::CoordType& max_distance) + { + return find_nearby_points(kdtree, center, max_distance, [](size_t) { + return true; + }); +} + + } // namespace Slic3r #endif /* slic3r_KDTreeIndirect_hpp_ */ diff --git a/src/libslic3r/Line.hpp b/src/libslic3r/Line.hpp index 751e59458..8631ee08b 100644 --- a/src/libslic3r/Line.hpp +++ b/src/libslic3r/Line.hpp @@ -113,7 +113,6 @@ public: Vector vector() const { return this->b - this->a; } Vector normal() const { return Vector((this->b(1) - this->a(1)), -(this->b(0) - this->a(0))); } bool intersection(const Line& line, Point* intersection) const; - double ccw(const Point& point) const { return point.ccw(*this); } // Clip a line with a bounding box. Returns false if the line is completely outside of the bounding box. bool clip_with_bbox(const BoundingBox &bbox); // Extend the line from both sides by an offset. diff --git a/src/libslic3r/Point.cpp b/src/libslic3r/Point.cpp index b2427d46d..9f56f96d6 100644 --- a/src/libslic3r/Point.cpp +++ b/src/libslic3r/Point.cpp @@ -108,36 +108,6 @@ bool Point::nearest_point(const Points &points, Point* point) const return true; } -/* Three points are a counter-clockwise turn if ccw > 0, clockwise if - * ccw < 0, and collinear if ccw = 0 because ccw is a determinant that - * gives the signed area of the triangle formed by p1, p2 and this point. - * In other words it is the 2D cross product of p1-p2 and p1-this, i.e. - * z-component of their 3D cross product. - * We return double because it must be big enough to hold 2*max(|coordinate|)^2 - */ -double Point::ccw(const Point &p1, const Point &p2) const -{ - static_assert(sizeof(coord_t) == 4, "Point::ccw() requires a 32 bit coord_t"); - return cross2((p2 - p1).cast(), (*this - p1).cast()); -// return cross2((p2 - p1).cast(), (*this - p1).cast()); -} - -double Point::ccw(const Line &line) const -{ - return this->ccw(line.a, line.b); -} - -// returns the CCW angle between this-p1 and this-p2 -// i.e. this assumes a CCW rotation from p1 to p2 around this -double Point::ccw_angle(const Point &p1, const Point &p2) const -{ - //FIXME this calculates an atan2 twice! Project one vector into the other! - double angle = atan2(p1.x() - (*this).x(), p1.y() - (*this).y()) - - atan2(p2.x() - (*this).x(), p2.y() - (*this).y()); - // we only want to return only positive angles - return angle <= 0 ? angle + 2*PI : angle; -} - Point Point::projection_onto(const MultiPoint &poly) const { Point running_projection = poly.first_point(); diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp index 6b430c2fe..84152ee9c 100644 --- a/src/libslic3r/Point.hpp +++ b/src/libslic3r/Point.hpp @@ -28,6 +28,9 @@ using Mat = Eigen::Matrix; template using Vec = Mat; +template +using DynVec = Eigen::Matrix; + // Eigen types, to replace the Slic3r's own types in the future. // Vector types with a fixed point coordinate base type. using Vec2crd = Eigen::Matrix; @@ -76,24 +79,35 @@ inline const auto &identity3d = identity<3, double>; inline bool operator<(const Vec2d &lhs, const Vec2d &rhs) { return lhs.x() < rhs.x() || (lhs.x() == rhs.x() && lhs.y() < rhs.y()); } -template -int32_t cross2(const Eigen::MatrixBase> &v1, const Eigen::MatrixBase> &v2) = delete; - -template -inline T cross2(const Eigen::MatrixBase> &v1, const Eigen::MatrixBase> &v2) -{ - return v1.x() * v2.y() - v1.y() * v2.x(); -} - +// Cross product of two 2D vectors. +// None of the vectors may be of int32_t type as the result would overflow. template inline typename Derived::Scalar cross2(const Eigen::MatrixBase &v1, const Eigen::MatrixBase &v2) { + static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector"); + static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector"); + static_assert(! std::is_same::value, "cross2(): Scalar type must not be int32_t, otherwise the cross product would overflow."); static_assert(std::is_same::value, "cross2(): Scalar types of 1st and 2nd operand must be equal."); return v1.x() * v2.y() - v1.y() * v2.x(); } -template -inline Eigen::Matrix perp(const Eigen::MatrixBase> &v) { return Eigen::Matrix(- v.y(), v.x()); } +// 2D vector perpendicular to the argument. +template +inline Eigen::Matrix perp(const Eigen::MatrixBase &v) +{ + static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "perp(): parameter is not a 2D vector"); + return { - v.y(), v.x() }; +} + +// Angle from v1 to v2, returning double atan2(y, x) normalized to <-PI, PI>. +template +inline double angle(const Eigen::MatrixBase &v1, const Eigen::MatrixBase &v2) { + static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "angle(): first parameter is not a 2D vector"); + static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "angle(): second parameter is not a 2D vector"); + auto v1d = v1.template cast(); + auto v2d = v2.template cast(); + return atan2(cross2(v1d, v2d), v1d.dot(v2d)); +} template Eigen::Matrix to_2d(const Eigen::MatrixBase> &ptN) { return { ptN.x(), ptN.y() }; } @@ -165,9 +179,6 @@ public: int nearest_point_index(const PointConstPtrs &points) const; int nearest_point_index(const PointPtrs &points) const; bool nearest_point(const Points &points, Point* point) const; - double ccw(const Point &p1, const Point &p2) const; - double ccw(const Line &line) const; - double ccw_angle(const Point &p1, const Point &p2) const; Point projection_onto(const MultiPoint &poly) const; Point projection_onto(const Line &line) const; }; diff --git a/src/libslic3r/Polygon.cpp b/src/libslic3r/Polygon.cpp index c7bbe9706..09f1c393d 100644 --- a/src/libslic3r/Polygon.cpp +++ b/src/libslic3r/Polygon.cpp @@ -171,50 +171,56 @@ Point Polygon::centroid() const return Point(Vec2d(c / (3. * area_sum))); } -// find all concave vertices (i.e. having an internal angle greater than the supplied angle) -// (external = right side, thus we consider ccw orientation) -Points Polygon::concave_points(double angle) const +// Filter points from poly to the output with the help of FilterFn. +// filter function receives two vectors: +// v1: this_point - previous_point +// v2: next_point - this_point +// and returns true if the point is to be copied to the output. +template +Points filter_points_by_vectors(const Points &poly, FilterFn filter) { - Points points; - angle = 2. * PI - angle + EPSILON; - - // check whether first point forms a concave angle - if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) <= angle) - points.push_back(this->points.front()); - - // check whether points 1..(n-1) form concave angles - for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++ p) - if (p->ccw_angle(*(p-1), *(p+1)) <= angle) - points.push_back(*p); - - // check whether last point forms a concave angle - if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) <= angle) - points.push_back(this->points.back()); - - return points; -} + // Last point is the first point visited. + Point p1 = poly.back(); + // Previous vector to p1. + Vec2d v1 = (p1 - *(poly.end() - 2)).cast(); -// find all convex vertices (i.e. having an internal angle smaller than the supplied angle) -// (external = right side, thus we consider ccw orientation) -Points Polygon::convex_points(double angle) const -{ - Points points; - angle = 2*PI - angle - EPSILON; - - // check whether first point forms a convex angle - if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) >= angle) - points.push_back(this->points.front()); - - // check whether points 1..(n-1) form convex angles - for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++p) { - if (p->ccw_angle(*(p-1), *(p+1)) >= angle) points.push_back(*p); + Points out; + for (Point p2 : poly) { + // p2 is next point to the currently visited point p1. + Vec2d v2 = (p2 - p1).cast(); + if (filter(v1, v2)) + out.emplace_back(p2); + v1 = v2; + p1 = p2; } - // check whether last point forms a convex angle - if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) >= angle) - points.push_back(this->points.back()); - - return points; + return out; +} + +template +Points filter_convex_concave_points_by_angle_threshold(const Points &poly, double angle_threshold, ConvexConcaveFilterFn convex_concave_filter) +{ + assert(angle_threshold >= 0.); + if (angle_threshold < EPSILON) { + double cos_angle = cos(angle_threshold); + return filter_points_by_vectors(poly, [convex_concave_filter, cos_angle](const Vec2d &v1, const Vec2d &v2){ + return convex_concave_filter(v1, v2) && v1.normalized().dot(v2.normalized()) < cos_angle; + }); + } else { + return filter_points_by_vectors(poly, [convex_concave_filter](const Vec2d &v1, const Vec2d &v2){ + return convex_concave_filter(v1, v2); + }); + } +} + +Points Polygon::convex_points(double angle_threshold) const +{ + return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) > 0.; }); +} + +Points Polygon::concave_points(double angle_threshold) const +{ + return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) < 0.; }); } // Projection of a point onto the polygon. diff --git a/src/libslic3r/Polygon.hpp b/src/libslic3r/Polygon.hpp index 089820565..77f3020be 100644 --- a/src/libslic3r/Polygon.hpp +++ b/src/libslic3r/Polygon.hpp @@ -65,8 +65,11 @@ public: void densify(float min_length, std::vector* lengths = nullptr); void triangulate_convex(Polygons* polygons) const; Point centroid() const; - Points concave_points(double angle = PI) const; - Points convex_points(double angle = PI) const; + // Considering CCW orientation of this polygon, find all convex resp. concave points + // with the angle at the vertex larger than a threshold. + // Zero angle_threshold means to accept all convex resp. concave points. + Points convex_points(double angle_threshold = 0.) const; + Points concave_points(double angle_threshold = 0.) const; // Projection of a point onto the polygon. Point point_projection(const Point &point) const; std::vector parameter_by_length() const; diff --git a/src/libslic3r/SLA/bicubic.h b/src/libslic3r/SLA/bicubic.h deleted file mode 100644 index 870d00dbd..000000000 --- a/src/libslic3r/SLA/bicubic.h +++ /dev/null @@ -1,186 +0,0 @@ -#ifndef BICUBIC_HPP -#define BICUBIC_HPP - -#include -#include -#include - -#include - -namespace Slic3r { - -namespace BicubicInternal { - // Linear kernel, to be able to test cubic methods with hat kernels. - template - struct LinearKernel - { - typedef T FloatType; - - static T a00() { return T(0.); } - static T a01() { return T(0.); } - static T a02() { return T(0.); } - static T a03() { return T(0.); } - static T a10() { return T(1.); } - static T a11() { return T(-1.); } - static T a12() { return T(0.); } - static T a13() { return T(0.); } - static T a20() { return T(0.); } - static T a21() { return T(1.); } - static T a22() { return T(0.); } - static T a23() { return T(0.); } - static T a30() { return T(0.); } - static T a31() { return T(0.); } - static T a32() { return T(0.); } - static T a33() { return T(0.); } - }; - - // Interpolation kernel aka Catmul-Rom aka Keyes kernel. - template - struct CubicCatmulRomKernel - { - typedef T FloatType; - - static T a00() { return 0; } - static T a01() { return (T)-0.5; } - static T a02() { return (T) 1.; } - static T a03() { return (T)-0.5; } - static T a10() { return (T) 1.; } - static T a11() { return 0; } - static T a12() { return (T)-5./2.; } - static T a13() { return (T) 3./2.; } - static T a20() { return 0; } - static T a21() { return (T) 0.5; } - static T a22() { return (T) 2.; } - static T a23() { return (T)-3./2.; } - static T a30() { return 0; } - static T a31() { return 0; } - static T a32() { return (T)-0.5; } - static T a33() { return (T) 0.5; } - }; - - // B-spline kernel - template - struct CubicBSplineKernel - { - typedef T FloatType; - - static T a00() { return (T) 1./6.; } - static T a01() { return (T) -3./6.; } - static T a02() { return (T) 3./6.; } - static T a03() { return (T) -1./6.; } - static T a10() { return (T) 4./6.; } - static T a11() { return 0; } - static T a12() { return (T) -6./6.; } - static T a13() { return (T) 3./6.; } - static T a20() { return (T) 1./6.; } - static T a21() { return (T) 3./6.; } - static T a22() { return (T) 3./6.; } - static T a23() { return (T)- 3./6.; } - static T a30() { return 0; } - static T a31() { return 0; } - static T a32() { return 0; } - static T a33() { return (T) 1./6.; } - }; - - template - inline T clamp(T a, T lower, T upper) - { - return (a < lower) ? lower : - (a > upper) ? upper : a; - } -} - -template -struct CubicKernel -{ - typedef typename KERNEL KernelInternal; - typedef typename KERNEL::FloatType FloatType; - - static FloatType kernel(FloatType x) - { - x = fabs(x); - if (x >= (FloatType)2.) - return 0.0f; - if (x <= (FloatType)1.) { - FloatType x2 = x * x; - FloatType x3 = x2 * x; - return KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3; - } - assert(x > (FloatType)1. && x < (FloatType)2.); - x -= (FloatType)1.; - FloatType x2 = x * x; - FloatType x3 = x2 * x; - return KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3; - } - - static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x) - { - const FloatType x2 = x*x; - const FloatType x3 = x*x*x; - return f0*(KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3) + - f1*(KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3) + - f2*(KERNEL::a20() + KERNEL::a21() * x + KERNEL::a22() * x2 + KERNEL::a23() * x3) + - f3*(KERNEL::a30() + KERNEL::a31() * x + KERNEL::a32() * x2 + KERNEL::a33() * x3); - } -}; - -// Linear splines -typedef CubicKernel> LinearKernelf; -typedef CubicKernel> LinearKerneld; -// Catmul-Rom splines -typedef CubicKernel> CubicCatmulRomKernelf; -typedef CubicKernel> CubicCatmulRomKerneld; -typedef CubicKernel> CubicInterpolationKernelf; -typedef CubicKernel> CubicInterpolationKerneld; -// Cubic B-splines -typedef CubicKernel> CubicBSplineKernelf; -typedef CubicKernel> CubicBSplineKerneld; - -template -static float cubic_interpolate(const Eigen::ArrayBase &F, const typename KERNEL::FloatType pt, const typename KERNEL::FloatType dx) -{ - typedef typename KERNEL::FloatType T; - const int w = int(F.size()); - const int ix = (int)floor(pt); - const T s = pt - (T)ix; - - if (ix > 1 && ix + 2 < w) { - // Inside the fully interpolated region. - return KERNEL::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s); - } - // Transition region. Extend with a constant function. - auto f = [&F, w](x) { return F[BicubicInternal::clamp(x, 0, w - 1)]; } - return KERNEL::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s); -} - -template -static float bicubic_interpolate(const Eigen::MatrixBase &F, const Eigen::Matrix &pt, const typename KERNEL::FloatType dx) -{ - typedef typename KERNEL::FloatType T; - const int w = F.cols(); - const int h = F.rows(); - const int ix = (int)floor(pt[0]); - const int iy = (int)floor(pt[1]); - const T s = pt[0] - (T)ix; - const T t = pt[1] - (T)iy; - - if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) { - // Inside the fully interpolated region. - return KERNEL::interpolate( - KERNEL::interpolate(F(ix-1,iy-1),F(ix ,iy-1),F(ix+1,iy-1),F(ix+2,iy-1),s), - KERNEL::interpolate(F(ix-1,iy ),F(ix ,iy ),F(ix+1,iy ),F(ix+2,iy ),s), - KERNEL::interpolate(F(ix-1,iy+1),F(ix ,iy+1),F(ix+1,iy+1),F(ix+2,iy+1),s), - KERNEL::interpolate(F(ix-1,iy+2),F(ix ,iy+2),F(ix+1,iy+2),F(ix+2,iy+2),s),t); - } - // Transition region. Extend with a constant function. - auto f = [&f, w, h](int x, int y) { return F(BicubicInternal::clamp(x,0,w-1),BicubicInternal::clamp(y,0,h-1)); } - return KERNEL::interpolate( - KERNEL::interpolate(f(ix-1,iy-1),f(ix ,iy-1),f(ix+1,iy-1),f(ix+2,iy-1),s), - KERNEL::interpolate(f(ix-1,iy ),f(ix ,iy ),f(ix+1,iy ),f(ix+2,iy ),s), - KERNEL::interpolate(f(ix-1,iy+1),f(ix ,iy+1),f(ix+1,iy+1),f(ix+2,iy+1),s), - KERNEL::interpolate(f(ix-1,iy+2),f(ix ,iy+2),f(ix+1,iy+2),f(ix+2,iy+2),s),t); -} - -} // namespace Slic3r - -#endif /* BICUBIC_HPP */ diff --git a/src/libslic3r/Subdivide.cpp b/src/libslic3r/Subdivide.cpp new file mode 100644 index 000000000..31988a851 --- /dev/null +++ b/src/libslic3r/Subdivide.cpp @@ -0,0 +1,218 @@ +#include "Subdivide.hpp" +#include "Point.hpp" + +namespace Slic3r{ + +indexed_triangle_set its_subdivide( + const indexed_triangle_set &its, float max_length) +{ + // same order as key order in Edge Divides + struct VerticesSequence + { + size_t start_index; + bool positive_order; + VerticesSequence(size_t start_index, bool positive_order = true) + : start_index(start_index), positive_order(positive_order){} + }; + // vertex index small, big vertex index from key.first to key.second + using EdgeDivides = std::map, VerticesSequence>; + struct Edges + { + Vec3f data[3]; + Vec3f lengths; + Edges(const Vec3crd &indices, const std::vector &vertices) + : lengths(-1.f,-1.f,-1.f) + { + const Vec3f &v0 = vertices[indices[0]]; + const Vec3f &v1 = vertices[indices[1]]; + const Vec3f &v2 = vertices[indices[2]]; + data[0] = v0 - v1; + data[1] = v1 - v2; + data[2] = v2 - v0; + } + float abs_sum(const Vec3f &v) + { + return abs(v[0]) + abs(v[1]) + abs(v[2]); + } + bool is_dividable(const float& max_length) { + Vec3f sum(abs_sum(data[0]), abs_sum(data[1]), abs_sum(data[2])); + Vec3i biggest_index = (sum[0] > sum[1]) ? + ((sum[0] > sum[2]) ? + ((sum[2] > sum[1]) ? + Vec3i(0, 2, 1) : + Vec3i(0, 1, 2)) : + Vec3i(2, 0, 1)) : + ((sum[1] > sum[2]) ? + ((sum[2] > sum[0]) ? + Vec3i(1, 2, 0) : + Vec3i(1, 0, 2)) : + Vec3i(2, 1, 0)); + for (int i = 0; i < 3; i++) { + int index = biggest_index[i]; + if (sum[index] <= max_length) return false; + lengths[index] = data[index].norm(); + if (lengths[index] <= max_length) continue; + + // calculate rest of lengths + for (int j = i + 1; j < 3; j++) { + index = biggest_index[j]; + lengths[index] = data[index].norm(); + } + return true; + } + return false; + } + }; + struct TriangleLengths + { + Vec3crd indices; + Vec3f l; // lengths + TriangleLengths(const Vec3crd &indices, const Vec3f &lengths) + : indices(indices), l(lengths) + {} + + int get_divide_index(float max_length) { + if (l[0] > l[1] && l[0] > l[2]) { + if (l[0] > max_length) return 0; + } else if (l[1] > l[2]) { + if (l[1] > max_length) return 1; + } else { + if (l[2] > max_length) return 2; + } + return -1; + } + + // divide triangle add new vertex to vertices + std::pair divide( + int divide_index, float max_length, + std::vector &vertices, + EdgeDivides &edge_divides) + { + // index to lengths and indices + size_t i0 = divide_index; + size_t i1 = (divide_index + 1) % 3; + size_t vi0 = indices[i0]; + size_t vi1 = indices[i1]; + std::pair key(vi0, vi1); + bool key_swap = false; + if (key.first > key.second) { + std::swap(key.first, key.second); + key_swap = true; + } + + float length = l[divide_index]; + size_t count_edge_vertices = static_cast(floor(length / max_length)); + float count_edge_segments = static_cast(count_edge_vertices + 1); + + auto it = edge_divides.find(key); + if (it == edge_divides.end()) { + // Create new vertices + VerticesSequence new_vs(vertices.size()); + Vec3f vf = vertices[key.first]; // copy + const Vec3f &vs = vertices[key.second]; + Vec3f dir = vs - vf; + for (size_t i = 1; i <= count_edge_vertices; ++i) { + float ratio = i / count_edge_segments; + vertices.push_back(vf + dir * ratio); + } + bool success; + std::tie(it,success) = edge_divides.insert({key, new_vs}); + assert(success); + } + const VerticesSequence &vs = it->second; + + int index_offset = count_edge_vertices/2; + size_t i2 = (divide_index + 2) % 3; + if (count_edge_vertices % 2 == 0 && key_swap == (l[i1] < l[i2])) { + --index_offset; + } + int sign = (vs.positive_order) ? 1 : -1; + size_t new_index = vs.start_index + sign*index_offset; + + size_t vi2 = indices[i2]; + const Vec3f &v2 = vertices[vi2]; + Vec3f new_edge = v2 - vertices[new_index]; + float new_len = new_edge.norm(); + + float ratio = (1 + index_offset) / count_edge_segments; + float len1 = l[i0] * ratio; + float len2 = l[i0] - len1; + if (key_swap) std::swap(len1, len2); + + Vec3crd indices1(vi0, new_index, vi2); + Vec3f lengths1(len1, new_len, l[i2]); + + Vec3crd indices2(new_index, vi1, vi2); + Vec3f lengths2(len2, l[i1], new_len); + + // append key for divided edge when neccesary + if (index_offset > 0) { + std::pair new_key(key.first, new_index); + bool new_key_swap = false; + if (new_key.first > new_key.second) { + std::swap(new_key.first, new_key.second); + new_key_swap = true; + } + if (edge_divides.find(new_key) == edge_divides.end()) { + // insert new + edge_divides.insert({new_key, (new_key_swap) ? + VerticesSequence(new_index - sign, !vs.positive_order) + : vs}); + } + } + + if (index_offset < int(count_edge_vertices)-1) { + std::pair new_key(new_index, key.second); + bool new_key_swap = false; + if (new_key.first > new_key.second) { + std::swap(new_key.first, new_key.second); + new_key_swap = true; + } + // bad order + if (edge_divides.find(new_key) == edge_divides.end()) { + edge_divides.insert({new_key, (new_key_swap) ? + VerticesSequence(vs.start_index + sign*(count_edge_vertices-1), !vs.positive_order) + : VerticesSequence(new_index + sign, vs.positive_order)}); + } + } + + return {TriangleLengths(indices1, lengths1), + TriangleLengths(indices2, lengths2)}; + } + }; + indexed_triangle_set result; + result.indices.reserve(its.indices.size()); + const std::vector &vertices = its.vertices; + result.vertices = vertices; // copy + std::queue tls; + + EdgeDivides edge_divides; + for (const Vec3crd &indices : its.indices) { + Edges edges(indices, vertices); + // speed up only sum not sqrt is apply + if (!edges.is_dividable(max_length)) { + // small triangle + result.indices.push_back(indices); + continue; + } + TriangleLengths tl(indices, edges.lengths); + do { + int divide_index = tl.get_divide_index(max_length); + if (divide_index < 0) { + // no dividing + result.indices.push_back(tl.indices); + if (tls.empty()) break; + tl = tls.front(); // copy + tls.pop(); + } else { + auto [tl1, tl2] = tl.divide(divide_index, max_length, + result.vertices, edge_divides); + tl = tl1; + tls.push(tl2); + } + } while (true); + } + return result; +} + +} diff --git a/src/libslic3r/Subdivide.hpp b/src/libslic3r/Subdivide.hpp new file mode 100644 index 000000000..f97e4b3c5 --- /dev/null +++ b/src/libslic3r/Subdivide.hpp @@ -0,0 +1,12 @@ +#ifndef libslic3r_Subdivide_hpp_ +#define libslic3r_Subdivide_hpp_ + +#include "TriangleMesh.hpp" + +namespace Slic3r { + +indexed_triangle_set its_subdivide(const indexed_triangle_set &its, float max_length); + +} + +#endif //libslic3r_Subdivide_hpp_ diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 5e0fb67b3..2131a9233 100644 --- a/src/libslic3r/libslic3r.h +++ b/src/libslic3r/libslic3r.h @@ -23,6 +23,13 @@ #include #include +#ifdef _WIN32 +// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation. +// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740 +// Thus it is recommended to use boost::container::deque instead. +#include +#endif // _WIN32 + #include "Technologies.hpp" #include "Semver.hpp" @@ -73,6 +80,16 @@ namespace Slic3r { extern Semver SEMVER; +// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation. +template> +using deque = +#ifdef _WIN32 + // Use boost implementation, which allocates blocks of 512 bytes instead of blocks of 8 bytes. + boost::container::deque; +#else // _WIN32 + std::deque; +#endif // _WIN32 + template inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); } diff --git a/src/libslic3r/pchheader.hpp b/src/libslic3r/pchheader.hpp index e6591f574..aeb79e3d8 100644 --- a/src/libslic3r/pchheader.hpp +++ b/src/libslic3r/pchheader.hpp @@ -64,6 +64,12 @@ #include #include #include +#ifdef _WIN32 +// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation. +// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740 +// Thus it is recommended to use boost::container::deque instead. +#include +#endif // _WIN32 #include #include #include diff --git a/t/geometry.t b/t/geometry.t index 0f37c0aa3..bb72b22e1 100644 --- a/t/geometry.t +++ b/t/geometry.t @@ -2,7 +2,7 @@ use Test::More; use strict; use warnings; -plan tests => 42; +plan tests => 30; BEGIN { use FindBin; @@ -189,63 +189,6 @@ my $polygons = [ is +Slic3r::Point->new(10, 0)->distance_to_line($line), 0, 'distance_to'; } -#========================================================== - -{ - my $square = Slic3r::Polygon->new_scale( - [100,100], - [200,100], - [200,200], - [100,200], - ); - is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in ccw square'; - is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in ccw square'; - - $square->make_clockwise; - is scalar(@{$square->concave_points(PI*4/3)}), 4, 'fuor concave vertices detected in cw square'; - is scalar(@{$square->convex_points(PI*2/3)}), 0, 'no convex vertices detected in cw square'; -} - -{ - my $square = Slic3r::Polygon->new_scale( - [150,100], - [200,100], - [200,200], - [100,200], - [100,100], - ); - is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon'; - is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square'; -} - -{ - my $square = Slic3r::Polygon->new_scale( - [200,200], - [100,200], - [100,100], - [150,100], - [200,100], - ); - is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon'; - is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square'; -} - -{ - my $triangle = Slic3r::Polygon->new( - [16000170,26257364], [714223,461012], [31286371,461008], - ); - is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle'; - is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle'; -} - -{ - my $triangle = Slic3r::Polygon->new( - [16000170,26257364], [714223,461012], [20000000,461012], [31286371,461012], - ); - is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle having collinear point'; - is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle having collinear point'; -} - { my $triangle = Slic3r::Polygon->new( [16000170,26257364], [714223,461012], [31286371,461008], diff --git a/t/perimeters.t b/t/perimeters.t index d3f96f122..adc2a6cec 100644 --- a/t/perimeters.t +++ b/t/perimeters.t @@ -223,7 +223,7 @@ use Slic3r::Test; if ($model eq 'cube_with_concave_hole') { # check that loop starts at a concave vertex - my $ccw_angle = $loop->first_point->ccw_angle(@$loop[-2,1]); + my $ccw_angle = $loop->[-2]->ccw($loop->first_point, $loop->[1]); my $convex = ($ccw_angle > PI); # whether the angle on the *right* side is convex $starts_on_convex_point = 1 if ($convex && $is_contour) || (!$convex && $is_hole); diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt index 69470aad1..248245182 100644 --- a/tests/libslic3r/CMakeLists.txt +++ b/tests/libslic3r/CMakeLists.txt @@ -8,6 +8,7 @@ add_executable(${_TEST_NAME}_tests test_clipper_utils.cpp test_color.cpp test_config.cpp + test_curve_fitting.cpp test_elephant_foot_compensation.cpp test_geometry.cpp test_placeholder_parser.cpp diff --git a/tests/libslic3r/test_curve_fitting.cpp b/tests/libslic3r/test_curve_fitting.cpp new file mode 100644 index 000000000..faf7839c7 --- /dev/null +++ b/tests/libslic3r/test_curve_fitting.cpp @@ -0,0 +1,118 @@ +#include +#include + +#include +#include +#include + +TEST_CASE("Curves: cubic b spline fit test", "[Curves]") { + using namespace Slic3r; + using namespace Slic3r::Geometry; + + auto fx = [&](size_t index) { + return float(index) / 200.0f; + }; + + auto fy = [&](size_t index) { + return 1.0f; + }; + + std::vector> observations { }; + std::vector observation_points { }; + std::vector weights { }; + for (size_t index = 0; index < 200; ++index) { + observations.push_back(Vec<1, float> { fy(index) }); + observation_points.push_back(fx(index)); + weights.push_back(1); + } + + Vec2f fmin { fx(0), fy(0) }; + Vec2f fmax { fx(200), fy(200) }; + + auto bspline = fit_cubic_bspline(observations, observation_points, weights, 1); + + Approx ap(1.0f); + ap.epsilon(0.1f); + + for (int p = 0; p < 200; ++p) { + float fitted_val = bspline.get_fitted_value(fx(p))(0); + float expected = fy(p); + + REQUIRE(fitted_val == ap(expected)); + + } +} + +TEST_CASE("Curves: quadratic f cubic b spline fit test", "[Curves]") { + using namespace Slic3r; + using namespace Slic3r::Geometry; + + auto fx = [&](size_t index) { + return float(index) / 100.0f; + }; + + auto fy = [&](size_t index) { + return (fx(index) - 1) * (fx(index) - 1); + }; + + std::vector> observations { }; + std::vector observation_points { }; + std::vector weights { }; + for (size_t index = 0; index < 200; ++index) { + observations.push_back(Vec<1, float> { fy(index) }); + observation_points.push_back(fx(index)); + weights.push_back(1); + } + + Vec2f fmin { fx(0), fy(0) }; + Vec2f fmax { fx(200), fy(200) }; + + auto bspline = fit_cubic_bspline(observations, observation_points, weights, 10); + + for (int p = 0; p < 200; ++p) { + float fitted_val = bspline.get_fitted_value(fx(p))(0); + float expected = fy(p); + + auto check = [](float a, float b) { + return abs(a - b) < 0.2f; + }; + //Note: checking is problematic, splines will not perfectly align + REQUIRE(check(fitted_val, expected)); + + } +} + +TEST_CASE("Curves: polynomial fit test", "[Curves]") { + using namespace Slic3r; + using namespace Slic3r::Geometry; + + auto fx = [&](size_t index) { + return float(index) / 100.0f; + }; + + auto fy = [&](size_t index) { + return (fx(index) - 1) * (fx(index) - 1); + }; + + std::vector> observations { }; + std::vector observation_points { }; + std::vector weights { }; + for (size_t index = 0; index < 200; ++index) { + observations.push_back(Vec<1, float> { fy(index) }); + observation_points.push_back(fx(index)); + weights.push_back(1); + } + + Vec2f fmin { fx(0), fy(0) }; + Vec2f fmax { fx(200), fy(200) }; + + Approx ap(1.0f); + ap.epsilon(0.1f); + + auto poly = fit_polynomial(observations, observation_points, weights, 2); + + REQUIRE(poly.coefficients(0, 0) == ap(1)); + REQUIRE(poly.coefficients(0, 1) == ap(-2)); + REQUIRE(poly.coefficients(0, 2) == ap(1)); +} + diff --git a/tests/libslic3r/test_geometry.cpp b/tests/libslic3r/test_geometry.cpp index 21f9a9663..34e625e6d 100644 --- a/tests/libslic3r/test_geometry.cpp +++ b/tests/libslic3r/test_geometry.cpp @@ -373,7 +373,43 @@ SCENARIO("Line distances", "[Geometry]"){ } } +SCENARIO("Calculating angles", "[Geometry]") +{ + GIVEN(("Vectors 30 degrees apart")) + { + std::vector> pts { + { {1000, 0}, { 866, 500 } }, + { { 866, 500 }, { 500, 866 } }, + { { 500, 866 }, { 0, 1000 } }, + { { -500, 866 }, { -866, 500 } } + }; + + THEN("Angle detected is 30 degrees") + { + for (auto &p : pts) + REQUIRE(is_approx(angle(p.first, p.second), M_PI / 6.)); + } + } + + GIVEN(("Vectors 30 degrees apart")) + { + std::vector> pts { + { { 866, 500 }, {1000, 0} }, + { { 500, 866 }, { 866, 500 } }, + { { 0, 1000 }, { 500, 866 } }, + { { -866, 500 }, { -500, 866 } } + }; + + THEN("Angle detected is -30 degrees") + { + for (auto &p : pts) + REQUIRE(is_approx(angle(p.first, p.second), - M_PI / 6.)); + } + } +} + SCENARIO("Polygon convex/concave detection", "[Geometry]"){ + static constexpr const double angle_threshold = M_PI / 3.; GIVEN(("A Square with dimension 100")){ auto square = Slic3r::Polygon /*new_scale*/(std::vector({ Point(100,100), @@ -381,13 +417,13 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(200,200), Point(100,200)})); THEN("It has 4 convex points counterclockwise"){ - REQUIRE(square.concave_points(PI*4/3).size() == 0); - REQUIRE(square.convex_points(PI*2/3).size() == 4); + REQUIRE(square.concave_points(angle_threshold).size() == 0); + REQUIRE(square.convex_points(angle_threshold).size() == 4); } THEN("It has 4 concave points clockwise"){ square.make_clockwise(); - REQUIRE(square.concave_points(PI*4/3).size() == 4); - REQUIRE(square.convex_points(PI*2/3).size() == 0); + REQUIRE(square.concave_points(angle_threshold).size() == 4); + REQUIRE(square.convex_points(angle_threshold).size() == 0); } } GIVEN("A Square with an extra colinearvertex"){ @@ -398,8 +434,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(100,200), Point(100,100)})); THEN("It has 4 convex points counterclockwise"){ - REQUIRE(square.concave_points(PI*4/3).size() == 0); - REQUIRE(square.convex_points(PI*2/3).size() == 4); + REQUIRE(square.concave_points(angle_threshold).size() == 0); + REQUIRE(square.convex_points(angle_threshold).size() == 4); } } GIVEN("A Square with an extra collinear vertex in different order"){ @@ -410,8 +446,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(150,100), Point(200,100)})); THEN("It has 4 convex points counterclockwise"){ - REQUIRE(square.concave_points(PI*4/3).size() == 0); - REQUIRE(square.convex_points(PI*2/3).size() == 4); + REQUIRE(square.concave_points(angle_threshold).size() == 0); + REQUIRE(square.convex_points(angle_threshold).size() == 4); } } @@ -422,8 +458,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(31286371,461008) })); THEN("it has three convex vertices"){ - REQUIRE(triangle.concave_points(PI*4/3).size() == 0); - REQUIRE(triangle.convex_points(PI*2/3).size() == 3); + REQUIRE(triangle.concave_points(angle_threshold).size() == 0); + REQUIRE(triangle.convex_points(angle_threshold).size() == 3); } } @@ -435,8 +471,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(31286371,461012) })); THEN("it has three convex vertices"){ - REQUIRE(triangle.concave_points(PI*4/3).size() == 0); - REQUIRE(triangle.convex_points(PI*2/3).size() == 3); + REQUIRE(triangle.concave_points(angle_threshold).size() == 0); + REQUIRE(triangle.convex_points(angle_threshold).size() == 3); } } GIVEN("A polygon with concave vertices with angles of specifically 4/3pi"){ @@ -453,8 +489,8 @@ SCENARIO("Polygon convex/concave detection", "[Geometry]"){ Point(38092663,692699),Point(52100125,692699) })); THEN("the correct number of points are detected"){ - REQUIRE(polygon.concave_points(PI*4/3).size() == 6); - REQUIRE(polygon.convex_points(PI*2/3).size() == 10); + REQUIRE(polygon.concave_points(angle_threshold).size() == 6); + REQUIRE(polygon.convex_points(angle_threshold).size() == 10); } } } diff --git a/xs/xsp/Line.xsp b/xs/xsp/Line.xsp index 777dc41fa..36181c3ba 100644 --- a/xs/xsp/Line.xsp +++ b/xs/xsp/Line.xsp @@ -41,7 +41,7 @@ Clone normal(); Clone vector(); double ccw(Point* point) - %code{% RETVAL = THIS->ccw(*point); %}; + %code{% RETVAL = cross2((THIS->a - *point).cast(), (THIS->b - THIS->a).cast()); %}; %{ Line* diff --git a/xs/xsp/Point.xsp b/xs/xsp/Point.xsp index beefc6249..2d70e0203 100644 --- a/xs/xsp/Point.xsp +++ b/xs/xsp/Point.xsp @@ -39,13 +39,7 @@ double perp_distance_to_line(Line* line) %code{% RETVAL = line->perp_distance_to(*THIS); %}; double ccw(Point* p1, Point* p2) - %code{% RETVAL = THIS->ccw(*p1, *p2); %}; - double ccw_angle(Point* p1, Point* p2) - %code{% RETVAL = THIS->ccw_angle(*p1, *p2); %}; - Point* projection_onto_polygon(Polygon* polygon) - %code{% RETVAL = new Point(THIS->projection_onto(*polygon)); %}; - Point* projection_onto_polyline(Polyline* polyline) - %code{% RETVAL = new Point(THIS->projection_onto(*polyline)); %}; + %code{% RETVAL = cross2((*p1 - *THIS).cast(), (*p2 - *p1).cast()); %}; Point* projection_onto_line(Line* line) %code{% RETVAL = new Point(THIS->projection_onto(*line)); %}; Point* negative() diff --git a/xs/xsp/Polygon.xsp b/xs/xsp/Polygon.xsp index a94425477..6b6d52524 100644 --- a/xs/xsp/Polygon.xsp +++ b/xs/xsp/Polygon.xsp @@ -39,8 +39,6 @@ %code{% THIS->triangulate_convex(&RETVAL); %}; Clone centroid(); Clone bounding_box(); - Points concave_points(double angle); - Points convex_points(double angle); Clone point_projection(Point* point) %code{% RETVAL = THIS->point_projection(*point); %}; Clone intersection(Line* line)