diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp index 217166f8c7..9b9c886ecc 100644 --- a/src/libslic3r/AABBTreeIndirect.hpp +++ b/src/libslic3r/AABBTreeIndirect.hpp @@ -446,6 +446,57 @@ namespace detail { } } + // Real-time collision detection, Ericson, Chapter 5 + template + static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c) + { + using Scalar = typename Vector::Scalar; + // Check if P in vertex region outside A + Vector ab = b - a; + Vector ac = c - a; + Vector ap = p - a; + Scalar d1 = ab.dot(ap); + Scalar d2 = ac.dot(ap); + if (d1 <= 0 && d2 <= 0) + return a; + // Check if P in vertex region outside B + Vector bp = p - b; + Scalar d3 = ab.dot(bp); + Scalar d4 = ac.dot(bp); + if (d3 >= 0 && d4 <= d3) + return b; + // Check if P in edge region of AB, if so return projection of P onto AB + Scalar vc = d1*d4 - d3*d2; + if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) { + Scalar v = d1 / (d1 - d3); + return a + v * ab; + } + // Check if P in vertex region outside C + Vector cp = p - c; + Scalar d5 = ab.dot(cp); + Scalar d6 = ac.dot(cp); + if (d6 >= 0 && d5 <= d6) + return c; + // Check if P in edge region of AC, if so return projection of P onto AC + Scalar vb = d5*d2 - d1*d6; + if (vb <= 0 && d2 >= 0 && d6 <= 0) { + Scalar w = d2 / (d2 - d6); + return a + w * ac; + } + // Check if P in edge region of BC, if so return projection of P onto BC + Scalar va = d3*d6 - d5*d4; + if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) { + Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); + return b + w * (c - b); + } + // P inside face region. Compute Q through its barycentric coordinates (u,v,w) + Scalar denom = Scalar(1.0) / (va + vb + vc); + Scalar v = vb * denom; + Scalar w = vc * denom; + return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w + }; + + // Nothing to do with COVID-19 social distancing. template struct IndexedTriangleSetDistancer { @@ -453,74 +504,36 @@ namespace detail { using IndexedFaceType = AIndexedFaceType; using TreeType = ATreeType; using VectorType = AVectorType; + using ScalarType = typename VectorType::Scalar; const std::vector &vertices; const std::vector &faces; const TreeType &tree; const VectorType origin; + + inline VectorType closest_point_to_origin(size_t primitive_index, + ScalarType& squared_distance){ + const auto &triangle = this->faces[primitive_index]; + VectorType closest_point = closest_point_to_triangle(origin, + this->vertices[triangle(0)].template cast(), + this->vertices[triangle(1)].template cast(), + this->vertices[triangle(2)].template cast()); + squared_distance = (origin - closest_point).squaredNorm(); + return closest_point; + } }; - // Real-time collision detection, Ericson, Chapter 5 - template - static inline Vector closest_point_to_triangle(const Vector &p, const Vector &a, const Vector &b, const Vector &c) - { - using Scalar = typename Vector::Scalar; - // Check if P in vertex region outside A - Vector ab = b - a; - Vector ac = c - a; - Vector ap = p - a; - Scalar d1 = ab.dot(ap); - Scalar d2 = ac.dot(ap); - if (d1 <= 0 && d2 <= 0) - return a; - // Check if P in vertex region outside B - Vector bp = p - b; - Scalar d3 = ab.dot(bp); - Scalar d4 = ac.dot(bp); - if (d3 >= 0 && d4 <= d3) - return b; - // Check if P in edge region of AB, if so return projection of P onto AB - Scalar vc = d1*d4 - d3*d2; - if (a != b && vc <= 0 && d1 >= 0 && d3 <= 0) { - Scalar v = d1 / (d1 - d3); - return a + v * ab; - } - // Check if P in vertex region outside C - Vector cp = p - c; - Scalar d5 = ab.dot(cp); - Scalar d6 = ac.dot(cp); - if (d6 >= 0 && d5 <= d6) - return c; - // Check if P in edge region of AC, if so return projection of P onto AC - Scalar vb = d5*d2 - d1*d6; - if (vb <= 0 && d2 >= 0 && d6 <= 0) { - Scalar w = d2 / (d2 - d6); - return a + w * ac; - } - // Check if P in edge region of BC, if so return projection of P onto BC - Scalar va = d3*d6 - d5*d4; - if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) { - Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); - return b + w * (c - b); - } - // P inside face region. Compute Q through its barycentric coordinates (u,v,w) - Scalar denom = Scalar(1.0) / (va + vb + vc); - Scalar v = vb * denom; - Scalar w = vc * denom; - return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w - }; - - template - static inline Scalar squared_distance_to_indexed_triangle_set_recursive( - IndexedTriangleSetDistancerType &distancer, + template + static inline Scalar squared_distance_to_indexed_primitives_recursive( + IndexedPrimitivesDistancerType &distancer, size_t node_idx, Scalar low_sqr_d, Scalar up_sqr_d, size_t &i, - Eigen::PlainObjectBase &c) + Eigen::PlainObjectBase &c) { - using Vector = typename IndexedTriangleSetDistancerType::VectorType; + using Vector = typename IndexedPrimitivesDistancerType::VectorType; if (low_sqr_d > up_sqr_d) return low_sqr_d; @@ -538,13 +551,9 @@ namespace detail { assert(node.is_valid()); if (node.is_leaf()) { - const auto &triangle = distancer.faces[node.idx]; - Vector c_candidate = closest_point_to_triangle( - distancer.origin, - distancer.vertices[triangle(0)].template cast(), - distancer.vertices[triangle(1)].template cast(), - distancer.vertices[triangle(2)].template cast()); - set_min((c_candidate - distancer.origin).squaredNorm(), node.idx, c_candidate); + Scalar sqr_dist; + Vector c_candidate = distancer.closest_point_to_origin(node.idx, sqr_dist); + set_min(sqr_dist, node.idx, c_candidate); } else { @@ -561,7 +570,7 @@ namespace detail { { size_t i_left; Vector c_left = c; - Scalar sqr_d_left = squared_distance_to_indexed_triangle_set_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left); + Scalar sqr_d_left = squared_distance_to_indexed_primitives_recursive(distancer, left_node_idx, low_sqr_d, up_sqr_d, i_left, c_left); set_min(sqr_d_left, i_left, c_left); looked_left = true; }; @@ -569,13 +578,13 @@ namespace detail { { size_t i_right; Vector c_right = c; - Scalar sqr_d_right = squared_distance_to_indexed_triangle_set_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right); + Scalar sqr_d_right = squared_distance_to_indexed_primitives_recursive(distancer, right_node_idx, low_sqr_d, up_sqr_d, i_right, c_right); set_min(sqr_d_right, i_right, c_right); looked_right = true; }; // must look left or right if in box - using BBoxScalar = typename IndexedTriangleSetDistancerType::TreeType::BoundingBox::Scalar; + using BBoxScalar = typename IndexedPrimitivesDistancerType::TreeType::BoundingBox::Scalar; if (node_left.bbox.contains(distancer.origin.template cast())) look_left(); if (node_right.bbox.contains(distancer.origin.template cast())) @@ -709,10 +718,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(); @@ -742,7 +756,7 @@ inline typename VectorType::Scalar squared_distance_to_indexed_triangle_set( auto distancer = detail::IndexedTriangleSetDistancer { vertices, faces, tree, point }; return tree.empty() ? Scalar(-1) : - detail::squared_distance_to_indexed_triangle_set_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits::infinity(), hit_idx_out, hit_point_out); + detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), std::numeric_limits::infinity(), hit_idx_out, hit_point_out); } // Decides if exists some triangle in defined radius on a 3D indexed triangle set using a pre-built AABBTreeIndirect::Tree. @@ -759,22 +773,22 @@ 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 { vertices, faces, tree, point }; - size_t hit_idx; - VectorType hit_point = VectorType::Ones() * (std::nan("")); + size_t hit_idx; + VectorType hit_point = VectorType::Ones() * (NaN); if(tree.empty()) { 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_primitives_recursive(distancer, size_t(0), Scalar(0), max_distance_squared, hit_idx, hit_point); return hit_point.allFinite(); } diff --git a/src/libslic3r/AABBTreeLines.hpp b/src/libslic3r/AABBTreeLines.hpp new file mode 100644 index 0000000000..7b9595419a --- /dev/null +++ b/src/libslic3r/AABBTreeLines.hpp @@ -0,0 +1,112 @@ +#ifndef SRC_LIBSLIC3R_AABBTREELINES_HPP_ +#define SRC_LIBSLIC3R_AABBTREELINES_HPP_ + +#include "libslic3r/Point.hpp" +#include "libslic3r/EdgeGrid.hpp" +#include "libslic3r/AABBTreeIndirect.hpp" +#include "libslic3r/Line.hpp" + +namespace Slic3r { + +namespace AABBTreeLines { + +namespace detail { + +template +struct IndexedLinesDistancer { + using LineType = ALineType; + using TreeType = ATreeType; + using VectorType = AVectorType; + using ScalarType = typename VectorType::Scalar; + + const std::vector &lines; + const TreeType &tree; + + const VectorType origin; + + inline VectorType closest_point_to_origin(size_t primitive_index, + ScalarType &squared_distance) { + VectorType nearest_point; + const LineType &line = lines[primitive_index]; + squared_distance = line_alg::distance_to_squared(line, origin, &nearest_point); + return nearest_point; + } +}; + +} + +// Build a balanced AABB Tree over a vector of float lines, balancing the tree +// on centroids of the lines. +// Epsilon is applied to the bounding boxes of the AABB Tree to cope with numeric inaccuracies +// during tree traversal. +template +inline AABBTreeIndirect::Tree<2, typename LineType::Scalar> build_aabb_tree_over_indexed_lines( + const std::vector &lines, + //FIXME do we want to apply an epsilon? + const float eps = 0) + { + using TreeType = AABBTreeIndirect::Tree<2, typename LineType::Scalar>; +// using CoordType = typename TreeType::CoordType; + using VectorType = typename TreeType::VectorType; + using BoundingBox = typename TreeType::BoundingBox; + + struct InputType { + size_t idx() const { + return m_idx; + } + const BoundingBox& bbox() const { + return m_bbox; + } + const VectorType& centroid() const { + return m_centroid; + } + + size_t m_idx; + BoundingBox m_bbox; + VectorType m_centroid; + }; + + std::vector input; + input.reserve(lines.size()); + const VectorType veps(eps, eps); + for (size_t i = 0; i < lines.size(); ++i) { + const LineType &line = lines[i]; + InputType n; + n.m_idx = i; + n.m_centroid = (line.a + line.b) * 0.5; + n.m_bbox = BoundingBox(line.a, line.a); + n.m_bbox.extend(line.b); + n.m_bbox.min() -= veps; + n.m_bbox.max() += veps; + input.emplace_back(n); + } + + TreeType out; + out.build(std::move(input)); + return out; +} + +// Finding a closest line, its closest point and squared distance to the closest point +// Returns squared distance to the closest point or -1 if the input is empty. +template +inline typename VectorType::Scalar squared_distance_to_indexed_lines( + const std::vector &lines, + const TreeType &tree, + const VectorType &point, + size_t &hit_idx_out, + Eigen::PlainObjectBase &hit_point_out) + { + using Scalar = typename VectorType::Scalar; + auto distancer = detail::IndexedLinesDistancer + { lines, tree, point }; + return tree.empty() ? + Scalar(-1) : + AABBTreeIndirect::detail::squared_distance_to_indexed_primitives_recursive(distancer, size_t(0), Scalar(0), + std::numeric_limits::infinity(), hit_idx_out, hit_point_out); +} + +} + +} + +#endif /* SRC_LIBSLIC3R_AABBTREELINES_HPP_ */ diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 2ebdb72e01..c6f01a65db 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -142,10 +142,12 @@ set(lisbslic3r_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 @@ -176,6 +178,8 @@ set(lisbslic3r_sources CustomGCode.hpp Arrange.hpp Arrange.cpp + NormalUtils.cpp + NormalUtils.hpp Orient.hpp Orient.cpp MultiPoint.cpp @@ -222,6 +226,8 @@ set(lisbslic3r_sources QuadricEdgeCollapse.cpp QuadricEdgeCollapse.hpp Semver.cpp + ShortEdgeCollapse.cpp + ShortEdgeCollapse.hpp ShortestPath.cpp ShortestPath.hpp SLAPrint.cpp @@ -264,6 +270,8 @@ set(lisbslic3r_sources Thread.hpp TriangleSelector.cpp TriangleSelector.hpp + TriangleSetSampling.cpp + TriangleSetSampling.hpp MTUtils.hpp VariableWidth.cpp VariableWidth.hpp diff --git a/src/libslic3r/ExtrusionEntity.cpp b/src/libslic3r/ExtrusionEntity.cpp index 2189f649fd..537adfdcc8 100644 --- a/src/libslic3r/ExtrusionEntity.cpp +++ b/src/libslic3r/ExtrusionEntity.cpp @@ -158,11 +158,10 @@ double ExtrusionLoop::length() const return len; } -bool ExtrusionLoop::split_at_vertex(const Point &point) +bool ExtrusionLoop::split_at_vertex(const Point &point, const double scaled_epsilon) { for (ExtrusionPaths::iterator path = this->paths.begin(); path != this->paths.end(); ++path) { - int idx = path->polyline.find_point(point); - if (idx != -1) { + if (int idx = path->polyline.find_point(point, scaled_epsilon); idx != -1) { if (this->paths.size() == 1) { // just change the order of points Polyline p1, p2; @@ -207,46 +206,57 @@ bool ExtrusionLoop::split_at_vertex(const Point &point) return false; } -std::pair ExtrusionLoop::get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const +ExtrusionLoop::ClosestPathPoint ExtrusionLoop::get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const { // Find the closest path and closest point belonging to that path. Avoid overhangs, if asked for. - size_t path_idx = 0; - Point p; - { - double min = std::numeric_limits::max(); - Point p_non_overhang; - size_t path_idx_non_overhang = 0; - double min_non_overhang = std::numeric_limits::max(); - for (const ExtrusionPath& path : this->paths) { - Point p_tmp = point.projection_onto(path.polyline); - double dist = (p_tmp - point).cast().norm(); - if (dist < min) { - p = p_tmp; - min = dist; - path_idx = &path - &this->paths.front(); - } - if (prefer_non_overhang && !is_bridge(path.role()) && dist < min_non_overhang) { - p_non_overhang = p_tmp; - min_non_overhang = dist; - path_idx_non_overhang = &path - &this->paths.front(); - } + ClosestPathPoint out{0, 0}; + double min2 = std::numeric_limits::max(); + ClosestPathPoint best_non_overhang{0, 0}; + double min2_non_overhang = std::numeric_limits::max(); + for (const ExtrusionPath &path : this->paths) { + std::pair foot_pt_ = foot_pt(path.polyline.points, point); + double d2 = (foot_pt_.second - point).cast().squaredNorm(); + if (d2 < min2) { + out.foot_pt = foot_pt_.second; + out.path_idx = &path - &this->paths.front(); + out.segment_idx = foot_pt_.first; + min2 = d2; } - if (prefer_non_overhang && min_non_overhang != std::numeric_limits::max()) { - // Only apply the non-overhang point if there is one. - path_idx = path_idx_non_overhang; - p = p_non_overhang; + if (prefer_non_overhang && !is_bridge(path.role()) && d2 < min2_non_overhang) { + best_non_overhang.foot_pt = foot_pt_.second; + best_non_overhang.path_idx = &path - &this->paths.front(); + best_non_overhang.segment_idx = foot_pt_.first; + min2_non_overhang = d2; } } - return std::make_pair(path_idx, p); + if (prefer_non_overhang && min2_non_overhang != std::numeric_limits::max()) + // Only apply the non-overhang point if there is one. + out = best_non_overhang; + return out; } // Splitting an extrusion loop, possibly made of multiple segments, some of the segments may be bridging. -void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang) +void ExtrusionLoop::split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon) { if (this->paths.empty()) return; - auto [path_idx, p] = get_closest_path_and_point(point, prefer_non_overhang); + auto [path_idx, segment_idx, p] = get_closest_path_and_point(point, prefer_non_overhang); + + // Snap p to start or end of segment_idx if closer than scaled_epsilon. + { + const Point *p1 = this->paths[path_idx].polyline.points.data() + segment_idx; + const Point *p2 = p1; + ++p2; + double d2_1 = (point - *p1).cast().squaredNorm(); + double d2_2 = (point - *p2).cast().squaredNorm(); + const double thr2 = scaled_epsilon * scaled_epsilon; + if (d2_1 < d2_2) { + if (d2_1 < thr2) p = *p1; + } else { + if (d2_2 < thr2) p = *p2; + } + } // now split path_idx in two parts const ExtrusionPath &path = this->paths[path_idx]; diff --git a/src/libslic3r/ExtrusionEntity.hpp b/src/libslic3r/ExtrusionEntity.hpp index 2fc6608000..6863a0df3b 100644 --- a/src/libslic3r/ExtrusionEntity.hpp +++ b/src/libslic3r/ExtrusionEntity.hpp @@ -314,9 +314,15 @@ public: const Point& last_point() const override { assert(this->first_point() == this->paths.back().polyline.points.back()); return this->first_point(); } Polygon polygon() const; double length() const override; - bool split_at_vertex(const Point &point); - void split_at(const Point &point, bool prefer_non_overhang); - std::pair get_closest_path_and_point(const Point& point, bool prefer_non_overhang) const; + bool split_at_vertex(const Point &point, const double scaled_epsilon = scaled(0.001)); + void split_at(const Point &point, bool prefer_non_overhang, const double scaled_epsilon = scaled(0.001)); + struct ClosestPathPoint + { + size_t path_idx; + size_t segment_idx; + Point foot_pt; + }; + ClosestPathPoint get_closest_path_and_point(const Point &point, bool prefer_non_overhang) const; void clip_end(double distance, ExtrusionPaths* paths) const; // Test, whether the point is extruded by a bridging flow. // This used to be used to avoid placing seams on overhangs, but now the EdgeGrid is used instead. diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index ec38f8e120..74b4035957 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -1567,7 +1567,8 @@ void GCode::_do_export(Print& print, GCodeOutputStream &file, ThumbnailsGenerato print.throw_if_canceled(); // Collect custom seam data from all objects. - m_seam_placer.init(print); + std::function throw_if_canceled_func = [&print]() { print.throw_if_canceled(); }; + m_seam_placer.init(print, throw_if_canceled_func); // BBS: priming logic is removed, always set first extruer here. //if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming)) @@ -2792,7 +2793,6 @@ GCode::LayerResult GCode::process_layer( m_wipe_tower->set_is_first_print(true); // 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) ? @@ -2981,9 +2981,9 @@ GCode::LayerResult GCode::process_layer( //This behaviour is same with cura if (is_infill_first && !first_layer) { 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 @@ -3161,14 +3161,11 @@ static std::unique_ptr calculate_layer_edge_grid(const Layer& la } -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, std::string 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); - //BBS: extrude contour of wall ccw, hole of wall cw, except spiral mode bool was_clockwise = loop.is_clockwise(); if (m_config.spiral_mode || !is_perimeter(loop.role())) @@ -3178,17 +3175,13 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou // 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_mode) { + if (!m_config.spiral_mode && description == "perimeter") { + assert(m_layer != nullptr); + bool is_outer_wall_first = m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill + || m_config.wall_infill_order == WallInfillOrder::InfillOuterInner; + m_seam_placer.place_seam(m_layer, loop, is_outer_wall_first, this->last_pos()); + } else loop.split_at(last_pos, false); - } - else { - //BBS - bool is_outer_wall_first = - m_config.wall_infill_order == WallInfillOrder::OuterInnerInfill || - m_config.wall_infill_order == WallInfillOrder::InfillOuterInner; - m_seam_placer.place_seam(loop, this->last_pos(), is_outer_wall_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 @@ -3298,14 +3291,14 @@ 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, std::string 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 ""; @@ -3327,24 +3320,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, "perimeter", -1.); } return gcode; } diff --git a/src/libslic3r/GCode.hpp b/src/libslic3r/GCode.hpp index 702f47fe73..29abf88def 100644 --- a/src/libslic3r/GCode.hpp +++ b/src/libslic3r/GCode.hpp @@ -325,8 +325,8 @@ private: std::string preamble(); // BBS std::string change_layer(coordf_t print_z, bool lazy_raise = false); - 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_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1.); + std::string extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1.); 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.); @@ -393,7 +393,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); diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index 97a64a64ef..bfb88d2f3c 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -1,987 +1,1384 @@ #include "SeamPlacer.hpp" +#include "tbb/parallel_for.h" +#include "tbb/blocked_range.h" +#include "tbb/parallel_reduce.h" +#include +#include +#include +#include + +#include "libslic3r/AABBTreeLines.hpp" +#include "libslic3r/KDTreeIndirect.hpp" #include "libslic3r/ExtrusionEntity.hpp" #include "libslic3r/Print.hpp" #include "libslic3r/BoundingBox.hpp" -#include "libslic3r/EdgeGrid.hpp" #include "libslic3r/ClipperUtils.hpp" -#include "libslic3r/SVG.hpp" #include "libslic3r/Layer.hpp" +#include "libslic3r/Geometry/Curves.hpp" +#include "libslic3r/ShortEdgeCollapse.hpp" +#include "libslic3r/TriangleSetSampling.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; - - - - -static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) +// ************ FOR BACKPORT COMPATIBILITY ONLY *************** +// Color mapping of a value into RGB false colors. +inline Vec3f value_to_rgbf(float minimum, float maximum, float value) { - // 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); - } + 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(); } +// *************************** +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) +// 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) { - 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; + 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); } - - -static Points::const_iterator project_point_to_polygon_and_insert(Polygon &polygon, const Point &pt, double eps) +float compute_angle_penalty(float ccw_angle) { - assert(polygon.points.size() >= 2); - if (polygon.points.size() <= 1) - if (polygon.points.size() == 1) - return polygon.points.begin(); + // 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)); +} - Point pt_min; - double d_min = std::numeric_limits::max(); - size_t i_min = size_t(-1); +/// Coordinate frame +class Frame +{ +public: + Frame() + { + mX = Vec3f(1, 0, 0); + mY = Vec3f(0, 1, 0); + mZ = Vec3f(0, 0, 1); + } - 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)); - } + Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) : mX(x), mY(y), mZ(z) {} + + 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); + } + + Vec3f to_world(const Vec3f &a) const { return a.x() * mX + a.y() * mY + a.z() * mZ; } + + 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, + const TriangleSetSamples & samples, + size_t negative_volumes_start_index) +{ + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: raycast visibility of " << samples.positions.size() << " samples over " << triangles.indices.size() << " triangles: end"; + + // prepare uniform samples of a hemisphere + float step_size = 1.0f / SeamPlacer::sqr_rays_per_sample_point; + std::vector precomputed_sample_directions(SeamPlacer::sqr_rays_per_sample_point * SeamPlacer::sqr_rays_per_sample_point); + for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_sample_point; ++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_sample_point; ++y_idx) { + size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_sample_point + 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(samples.positions.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, &samples](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 s_idx = r.begin(); s_idx < r.end(); ++s_idx) { + result[s_idx] = 1.0f; + constexpr float decrease_step = 1.0f / (SeamPlacer::sqr_rays_per_sample_point * SeamPlacer::sqr_rays_per_sample_point); + const Vec3f ¢er = samples.positions[s_idx]; + const Vec3f &normal = samples.normals[s_idx]; + // 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(); + 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.01f).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 && its_face_normal(triangles, hitpoint.id).dot(final_ray_dir) <= 0) { result[s_idx] -= decrease_step; } + } else { // TODO improve logic for order based boolean operations - consider order of volumes + bool casting_from_negative_volume = samples.triangle_indices[s_idx] >= negative_volumes_start_index; - // 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); - 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); - } - else if (! external_first) { - // Internal perimeter printed before the external. - // First get list of external seams. - std::vector ext_seams; - for (size_t i = 0; i < m_plan.size(); ++i) { - if (m_plan[i].external) - ext_seams.emplace_back(i); - } - - 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; + Vec3d ray_origin_d = (center + normal * 0.01f).cast(); // start above surface. + if (casting_from_negative_volume) { // 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.01f).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 counter = 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 + counter -= 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 { + counter += sgn(face_normal.dot(final_ray_dir)); } } + if (counter == 0) { result[s_idx] -= decrease_step; } } } - - // Only accept seam that is reasonably close. - double limit_dist_sqr = std::pow(double(scale_((ext_seam_idx - m_plan_idx) * nozzle_diameter * 2.)), 2.); - if (ext_seam_idx != size_t(-1) && 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]); - 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(); - } } } - 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; - } - } + }); + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: raycast visibility of " << samples.positions.size() << " samples over " << triangles.indices.size() << " triangles: end"; - // 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_idx; + return result; } -constexpr float CLOSE_TO_LAST_SEAM_THRESHOLD = 5; - -// Returns a seam for an EXTERNAL 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) +std::vector calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector &lengths, float min_arm_length) { - assert(loop.role() == erExternalPerimeter); - Polygon polygon = loop.polygon(); - bool was_clockwise = polygon.make_counter_clockwise(); - BoundingBox polygon_bb = polygon.bounding_box(); - const coord_t nozzle_r = coord_t(scale_(0.5 * nozzle_dmr) + 0.5); + std::vector result(polygon.size()); - size_t po_idx = std::find(m_po_list.begin(), m_po_list.end(), po) - m_po_list.begin(); + if (polygon.size() == 1) { result[0] = 0.0f; } - // 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->get_layer_at_printz(layer.print_z); - m_last_po = po; - m_last_print_z = layer.print_z; - m_last_layer_po = layer_po; - } - if (! layer_po) - return last_pos; + size_t idx_prev = 0; + size_t idx_curr = 0; + size_t idx_next = 0; - // Index of this layer in the respective PrintObject. - size_t layer_idx = layer_po->id() - po->layers().front()->id(); // raft layers + float distance_to_prev = 0; + float distance_to_next = 0; - assert(layer_idx < po->layer_count()); - - if (this->is_custom_seam_on_layer(layer_idx, po_idx)) { - // 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); + // 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]; } - if (seam_position != spRandom) { - // Retrieve the last start position for this object. - float last_pos_weight = 1.f; - - 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, polygon_bb); - if (pos.has_value()) { - last_pos = *pos; - last_pos_weight = is_custom_enforcer_on_layer(layer_idx, po_idx) ? 0.f : 1.f; - } - } - } - 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 - - // BBS: if the project point of current nozzle position is close to the last seam of this object - // then we think the current nozzle position is almost same with last same. - // So that last seam can be treat as one factor even in cost based strategy to make seam more posible to be aligned - if (po != nullptr) { - std::optional pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, polygon_bb); - if (pos.has_value()) { - Point project_point = polygon.point_projection(last_pos); - if ((pos->cast() - project_point.cast()).squaredNorm() < std::pow(scale_(CLOSE_TO_LAST_SEAM_THRESHOLD), 2.)) - last_pos = *pos; - } - } + 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()); } - // 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(); + // 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()); } - // Parametrize the polygon by its length. - std::vector lengths = polygon.parameter_by_length(); + // 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)); - // 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 = 5.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); - } + // increase idx_curr by one + float curr_distance = lengths[idx_curr]; + idx_curr++; + distance_to_prev += curr_distance; + distance_to_next -= curr_distance; + } - // 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. - 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; - } - } - - if (seam_position == spAligned || seam_position == spNearest) - m_seam_history.add_seam(po, polygon.points[idx_min], polygon_bb); - - - // 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); + return result; } - -Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_idx, - bool* saw_custom) const +struct CoordinateFunctor { - // Parametrize the polygon by its length. - const std::vector lengths = polygon.parameter_by_length(); + const std::vector *coordinates; + CoordinateFunctor(const std::vector *coords) : coordinates(coords) {} + CoordinateFunctor() : coordinates(nullptr) {} - // 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); + const float &operator()(size_t idx, size_t dim) const { return coordinates->operator[](idx)[dim]; } +}; - bool has_enforcers = ! enforcers_idxs.empty(); - bool has_blockers = ! blockers_idxs.empty(); - if (saw_custom) - *saw_custom = has_enforcers || has_blockers; +// structure to store global information about the model - occlusion hits, enforcers, blockers +struct GlobalModelInfo +{ + TriangleSetSamples mesh_samples; + std::vector mesh_samples_visibility; + CoordinateFunctor mesh_samples_coordinate_functor; + KDTreeIndirect<3, float, CoordinateFunctor> mesh_samples_tree{CoordinateFunctor{}}; + float mesh_samples_radius; - assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end())); - assert(std::is_sorted(blockers_idxs.begin(), blockers_idxs.end())); - std::vector edges; + indexed_triangle_set enforcers; + indexed_triangle_set blockers; + AABBTreeIndirect::Tree<3, float> enforcers_tree; + AABBTreeIndirect::Tree<3, float> blockers_tree; - // 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 + bool is_enforced(const Vec3f &position, float radius) const { - // 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; + 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 + { + std::vector points = find_nearby_points(mesh_samples_tree, position, mesh_samples_radius); + if (points.empty()) { return 1.0f; } + + auto compute_dist_to_plane = [](const Vec3f &position, const Vec3f &plane_origin, const Vec3f &plane_normal) { + Vec3f orig_to_point = position - plane_origin; + return std::abs(orig_to_point.dot(plane_normal)); + }; + + float total_weight = 0; + float total_visibility = 0; + for (size_t i = 0; i < points.size(); ++i) { + size_t sample_idx = points[i]; + + Vec3f sample_point = this->mesh_samples.positions[sample_idx]; + Vec3f sample_normal = this->mesh_samples.normals[sample_idx]; + + float weight = mesh_samples_radius - compute_dist_to_plane(position, sample_point, sample_normal); + weight += (mesh_samples_radius - (position - sample_point).norm()); + total_visibility += weight * mesh_samples_visibility[sample_idx]; + total_weight += weight; } - // 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 total_visibility / total_weight; + } + +#ifdef DEBUG_FILES + void debug_export(const indexed_triangle_set &obj_mesh) const + { + indexed_triangle_set divided_mesh = obj_mesh; + Slic3r::CNumericLocalesSetter locales_setter; + + { + auto filename = debug_out_path("visiblity.obj"); + FILE *fp = boost::nowide::fopen(filename.c_str(), "w"); + if (fp == nullptr) { + BOOST_LOG_TRIVIAL(error) << "stl_write_obj: Couldn't open " << filename << " 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); + } + + { + auto filename = debug_out_path("visiblity_samples.obj"); + FILE *fp = boost::nowide::fopen(filename.c_str(), "w"); + if (fp == nullptr) { + BOOST_LOG_TRIVIAL(error) << "stl_write_obj: Couldn't open " << filename << " for writing"; + return; + } + + for (size_t i = 0; i < mesh_samples.positions.size(); ++i) { + float visibility = mesh_samples_visibility[i]; + Vec3f color = value_to_rgbf(0.0f, 1.0f, visibility); + fprintf(fp, "v %f %f %f %f %f %f\n", mesh_samples.positions[i](0), mesh_samples.positions[i](1), mesh_samples.positions[i](2), color(0), color(1), color(2)); + } + 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) { + ExtrusionRole role = perimeter->role(); + if (perimeter->is_loop()) { + for (const ExtrusionPath &path : static_cast(perimeter)->paths) { + if (path.role() == ExtrusionRole::erExternalPerimeter) { role = ExtrusionRole::erExternalPerimeter; } + } + } + + if (role == ExtrusionRole::erExternalPerimeter || + (is_perimeter(role) && 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); } } - return running_total; + } + + 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); + } + + 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; } + Polygon polygon = orig_polygon; + bool was_clockwise = polygon.make_counter_clockwise(); + + std::vector lengths{}; + for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) { lengths.push_back((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm()); } + lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.1)); + std::vector polygon_angles = calculate_polygon_angles_at_vertices(polygon, lengths, SeamPlacer::polygon_local_angles_arm_distance); + + result.perimeters.push_back({}); + Perimeter &perimeter = result.perimeters.back(); + + 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 ? -polygon_angles[orig_angle_index] : polygon_angles[orig_angle_index]; + orig_angle_index++; + orig_point = true; + } + + if (global_model_info.is_enforced(position, perimeter.flow_width)) { type = EnforcedBlockedSeamPoint::Enforced; } + + if (global_model_info.is_blocked(position, perimeter.flow_width)) { type = EnforcedBlockedSeamPoint::Blocked; } + some_point_enforced = some_point_enforced || type == EnforcedBlockedSeamPoint::Enforced; + + 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; + } + } + } + + result.points.emplace_back(position, perimeter, local_ccw_angle, type); + } + + perimeter.end_index = result.points.size(); + + if (some_point_enforced) { + // We will patches of enforced points (patch: continuous section of enforced points), choose + // the longest patch, and select the middle point or sharp point (depending on the angle) + // this point will have high priority on this perimeter + size_t perimeter_size = perimeter.end_index - perimeter.start_index; + const auto next_index = [&](size_t idx) { return perimeter.start_index + Slic3r::next_idx_modulo(idx - perimeter.start_index, perimeter_size); }; + + std::vector patches_starts_ends; + for (size_t i = perimeter.start_index; i < perimeter.end_index; ++i) { + if (result.points[i].type != EnforcedBlockedSeamPoint::Enforced && result.points[next_index(i)].type == EnforcedBlockedSeamPoint::Enforced) { + patches_starts_ends.push_back(next_index(i)); + } + if (result.points[i].type == EnforcedBlockedSeamPoint::Enforced && result.points[next_index(i)].type != EnforcedBlockedSeamPoint::Enforced) { + patches_starts_ends.push_back(next_index(i)); + } + } + // if patches_starts_ends are empty, it means that the whole perimeter is enforced.. don't do anything in that case + if (!patches_starts_ends.empty()) { + // if the first point in the patches is not enforced, it marks a patch end. in that case, put it to the end and start on next + // to simplify the processing + assert(patches_starts_ends.size() % 2 == 0); + bool start_on_second = false; + if (result.points[patches_starts_ends[0]].type != EnforcedBlockedSeamPoint::Enforced) { + start_on_second = true; + patches_starts_ends.push_back(patches_starts_ends[0]); + } + // now pick the longest patch + std::pair longest_patch{0, 0}; + auto patch_len = [perimeter_size](const std::pair &start_end) { + if (start_end.second < start_end.first) { + return start_end.first + (perimeter_size - start_end.second); + } else { + return start_end.second - start_end.first; + } + }; + for (size_t patch_idx = start_on_second ? 1 : 0; patch_idx < patches_starts_ends.size(); patch_idx += 2) { + std::pair current_patch{patches_starts_ends[patch_idx], patches_starts_ends[patch_idx + 1]}; + if (patch_len(longest_patch) < patch_len(current_patch)) { longest_patch = current_patch; } + } + std::vector viable_points_indices; + std::vector large_angle_points_indices; + for (size_t point_idx = longest_patch.first; point_idx != longest_patch.second; point_idx = next_index(point_idx)) { + viable_points_indices.push_back(point_idx); + if (std::abs(result.points[point_idx].local_ccw_angle) > SeamPlacer::sharp_angle_snapping_threshold) { large_angle_points_indices.push_back(point_idx); } + } + assert(viable_points_indices.size() > 0); + if (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 { + size_t central_idx = large_angle_points_indices.size() / 2; + result.points[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 - 1) { + // 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, std::function throw_if_canceled) +{ + 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); + } + } + } + throw_if_canceled(); + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: gather occlusion meshes: end"; + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: decimate: start"; + its_short_edge_collpase(triangle_set, 25000); + its_short_edge_collpase(negative_volumes_set, 25000); + + 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: decimate: end"; + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: Compute visibility sample points: start"; + + result.mesh_samples = sample_its_uniform_parallel(SeamPlacer::raycasting_visibility_samples_count, triangle_set); + result.mesh_samples_coordinate_functor = CoordinateFunctor(&result.mesh_samples.positions); + result.mesh_samples_tree = KDTreeIndirect<3, float, CoordinateFunctor>(result.mesh_samples_coordinate_functor, result.mesh_samples.positions.size()); + + // The following code determines search area for random visibility samples on the mesh when calculating visibility of each perimeter point + // number of random samples in the given radius (area) is approximately poisson distribution + // to compute ideal search radius (area), we use exponential distribution (complementary distr to poisson) + // parameters of exponential distribution to compute area that will have with probability="probability" more than given number of samples="samples" + float probability = 0.9f; + float samples = 4; + float density = SeamPlacer::raycasting_visibility_samples_count / result.mesh_samples.total_area; + // exponential probability distrubtion function is : f(x) = P(X > x) = e^(l*x) where l is the rate parameter (computed as 1/u where u is mean value) + // probability that sampled area A with S samples contains more than samples count: + // P(S > samples in A) = e^-(samples/(density*A)); express A: + float search_area = samples / (-logf(probability) * density); + float search_radius = sqrt(search_area / PI); + result.mesh_samples_radius = search_radius; + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: Compute visiblity sample points: end"; + throw_if_canceled(); + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: Mesh sample raidus: " << result.mesh_samples_radius; + + 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); + + throw_if_canceled(); + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: build AABB tree: end"; + result.mesh_samples_visibility = raycast_visibility(raycasting_tree, triangle_set, result.mesh_samples, negative_volumes_start_index); + throw_if_canceled(); +#ifdef DEBUG_FILES + result.debug_export(triangle_set); +#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(); + + for (const ModelVolume *mv : po->model_object()->volumes) { + if (mv->is_seam_painted()) { + auto model_transformation = obj_transform * mv->get_matrix(); + + indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER); + its_transform(enforcers, model_transformation); + its_merge(result.enforcers, enforcers); + + indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER); + its_transform(blockers, model_transformation); + its_merge(result.blockers, blockers); + } + } + + 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); + + BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end"; +} + +struct SeamComparator +{ + SeamPosition setup; + float angle_importance; + explicit SeamComparator(SeamPosition setup) : setup(setup) + { + angle_importance = setup == spNearest ? SeamPlacer::angle_importance_nearest : SeamPlacer::angle_importance_aligned; + } + + // 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; } + + // Blockers/Enforcers discrimination, top priority + 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 0.5 mm inside) + if (a.embedded_distance < -0.5f && b.embedded_distance > -0.5f) { return true; } + if (b.embedded_distance < -0.5f && a.embedded_distance > -0.5f) { return false; } + + if (setup == SeamPosition::spRear && a.position.y() != b.position.y()) { return a.position.y() > b.position.y(); } + + float distance_penalty_a = 0.0f; + float distance_penalty_b = 0.0f; + if (setup == spNearest) { + distance_penalty_a = 1.0f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + distance_penalty_b = 1.0f - gauss((b.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + } + + // the penalites are kept close to range [0-1.x] however, it should not be relied upon + float penalty_a = a.overhang + a.visibility + angle_importance * compute_angle_penalty(a.local_ccw_angle) + distance_penalty_a; + float penalty_b = b.overhang + b.visibility + 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) && abs(a.overhang - b.overhang) > (0.1f * a.perimeter.flow_width)) { return a.overhang < b.overhang; } + + // prefer hidden points (more than 0.5 mm inside) + if (a.embedded_distance < -0.5f && b.embedded_distance > -0.5f) { return true; } + if (b.embedded_distance < -0.5f && a.embedded_distance > -0.5f) { return false; } + + if (setup == SeamPosition::spRandom) { return true; } + + if (setup == SeamPosition::spRear) { return a.position.y() + SeamPlacer::seam_align_score_tolerance * 5.0f > b.position.y(); } + + float penalty_a = a.overhang + a.visibility + angle_importance * compute_angle_penalty(a.local_ccw_angle); + float penalty_b = b.overhang + b.visibility + 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); } +}; + +#ifdef DEBUG_FILES +void debug_export_points(const std::vector &layers, const BoundingBox &bounding_box, const SeamComparator &comparator) +{ + for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) { + std::string angles_file_name = debug_out_path(("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, -compute_angle_penalty(point.local_ccw_angle)); + max_weight = std::max(max_weight, -compute_angle_penalty(point.local_ccw_angle)); + } + + std::string visiblity_file_name = debug_out_path(("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(("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(("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, -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 - 1 ? 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 - 1 ? 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; +} + +class PerimeterDistancer +{ + std::vector lines; + AABBTreeIndirect::Tree<2, double> tree; + +public: + PerimeterDistancer(const Layer *layer) + { + static const float eps = float(scale_(g_config_slice_closing_radius)); + // merge with offset + ExPolygons merged = layer->merged(eps); + // ofsset back + ExPolygons layer_outline = offset_ex(merged, -eps); + for (const ExPolygon &island : layer_outline) { + assert(island.contour.is_counter_clockwise()); + for (const auto &line : island.contour.lines()) { lines.emplace_back(unscale(line.a), unscale(line.b)); } + for (const Polygon &hole : island.holes) { + assert(hole.is_clockwise()); + for (const auto &line : hole.lines()) { lines.emplace_back(unscale(line.a), unscale(line.b)); } + } + } + tree = AABBTreeLines::build_aabb_tree_over_indexed_lines(lines); + } + + float distance_from_perimeter(const Point &point) const + { + Vec2d p = unscale(point); + size_t hit_idx_out{}; + Vec2d hit_point_out = Vec2d::Zero(); + auto distance = AABBTreeLines::squared_distance_to_indexed_lines(lines, tree, p, hit_idx_out, hit_point_out); + if (distance < 0) { return std::numeric_limits::max(); } + + distance = sqrt(distance); + const Linef &line = lines[hit_idx_out]; + Vec2d v1 = line.b - line.a; + Vec2d v2 = p - line.a; + if ((v1.x() * v2.y()) - (v1.y() * v2.x()) > 0.0) { distance *= -1; } + return distance; + } +}; + +} // 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_distancer; + if (r.begin() > 0) { // previous layer exists + prev_layer_distancer = std::make_unique(po->layers()[r.begin() - 1]); + } + + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + size_t regions_with_perimeter = 0; + for (const LayerRegion *region : po->layers()[layer_idx]->regions()) { + if (region->perimeters.entities.size() > 0) { regions_with_perimeter++; } + }; + bool should_compute_layer_embedding = regions_with_perimeter > 1; + std::unique_ptr current_layer_distancer = std::make_unique(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_distancer.get() != nullptr) { + perimeter_point.overhang = (prev_layer_distancer->distance_from_perimeter(point) + 0.5f * perimeter_point.perimeter.flow_width - + tan(SeamPlacer::overhang_angle_threshold) * po->layers()[layer_idx]->height) / + (3.0f * perimeter_point.perimeter.flow_width); + // NOTE disables the feature to place seams on slowly decreasing areas. Remove the following line to enable. + perimeter_point.overhang = perimeter_point.overhang < 0.0f ? 0.0f : perimeter_point.overhang; + } + + if (should_compute_layer_embedding) { // search for embedded perimeter points (points hidden inside the print ,e.g. multimaterial join, best position for seam) + perimeter_point.embedded_distance = current_layer_distancer->distance_from_perimeter(point) + 0.5f * perimeter_point.perimeter.flow_width; + } + } + + prev_layer_distancer.swap(current_layer_distancer); + } + }); +} + +// 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(). +std::optional> SeamPlacer::find_next_seam_in_layer(const std::vector &layers, + const Vec3f & projected_position, + const size_t layer_idx, + const float max_distance, + const SeamPlacerImpl::SeamComparator & comparator) const +{ + using namespace SeamPlacerImpl; + std::vector nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position, max_distance); + + if (nearby_points_indices.empty()) { return {}; } + + 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 {}; + } + + // 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 * max_distance)) { + return {std::pair{layer_idx, nearest_point.perimeter.seam_index}}; + } + + // First try to align the nearest, then try the best nearby + if (comparator.is_first_not_much_worse(nearest_point, next_layer_seam)) { return {std::pair{layer_idx, nearest_point_index}}; } + // 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)) { return {std::pair{layer_idx, best_nearby_point_index}}; } + + return {}; +} + +std::vector> SeamPlacer::find_seam_string(const PrintObject * po, + std::pair start_seam, + const SeamPlacerImpl::SeamComparator &comparator) const +{ + const std::vector &layers = m_seam_per_object.find(po)->second.layers; + int layer_idx = start_seam.first; + + // initialize searching for seam string - cluster of nearby seams on previous and next layers + int next_layer = layer_idx + 1; + int step = 1; + std::pair prev_point_index = start_seam; + std::vector> seam_string{start_seam}; + + auto reverse_lookup_direction = [&]() { + step = -1; + prev_point_index = start_seam; + next_layer = layer_idx - 1; }; - // 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); + while (next_layer >= 0) { + if (next_layer >= int(layers.size())) { + reverse_lookup_direction(); + if (next_layer < 0) { break; } + } + float max_distance = SeamPlacer::seam_align_tolerable_dist_factor * layers[start_seam.first].points[start_seam.second].perimeter.flow_width; + Vec3f prev_position = layers[prev_point_index.first].points[prev_point_index.second].position; + Vec3f projected_position = prev_position; + projected_position.z() = float(po->get_layer(next_layer)->slice_z); - 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(); - } - 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; + std::optional> maybe_next_seam = find_next_seam_in_layer(layers, projected_position, next_layer, max_distance, comparator); + + if (maybe_next_seam.has_value()) { + // For old macOS (pre 10.14), std::optional does not have .value() method, so the code is using operator*() instead. + seam_string.push_back(maybe_next_seam.operator*()); + prev_point_index = seam_string.back(); + // String added, prev_point_index updated + } else { + if (step == 1) { + reverse_lookup_direction(); + if (next_layer < 0) { break; } + } else { break; } } + next_layer += step; } - - 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]; - } + return seam_string; } - - - - - - - -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 +// 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) { - enforcers_idxs.clear(); - blockers_idxs.clear(); + using namespace SeamPlacerImpl; - 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 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 t_e) ? t_s + half_dist : t_e - half_dist; - - 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; -} - - - -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)) + // Prepares Debug files for writing. +#ifdef DEBUG_FILES + Slic3r::CNumericLocalesSetter locales_setter; + auto clusters_f = debug_out_path("seam_clusters.obj"); + 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; - - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs); - - 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; - } - } - -#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"); - - if (! blockers_idxs.empty()) { - ; - } - - - size_t min_idx = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); - - for (size_t i=0; i 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; + // 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; + } } - std::optional out; + // sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice + std::stable_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]); + }); - auto seams_it = m_data_last_layer.find(po); - if (seams_it == m_data_last_layer.end()) - return out; + // 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> alternative_seam_string; + std::vector observations; + std::vector observation_points; + std::vector weights; - 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. + int global_index = 0; + while (global_index < int(seams.size())) { + size_t layer_idx = seams[global_index].first; + size_t seam_index = seams[global_index].second; + global_index++; + 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 { + seam_string = this->find_seam_string(po, {layer_idx, seam_index}, comparator); + size_t step_size = 1 + seam_string.size() / 20; + for (size_t alternative_start = 0; alternative_start < seam_string.size(); alternative_start += step_size) { + size_t start_layer_idx = seam_string[alternative_start].first; + size_t seam_idx = layers[start_layer_idx].points[seam_string[alternative_start].second].perimeter.seam_index; + alternative_seam_string = this->find_seam_string(po, std::pair(start_layer_idx, seam_idx), comparator); + if (alternative_seam_string.size() > seam_string.size()) { seam_string = std::move(alternative_seam_string); } + } + if (seam_string.size() < seam_align_minimum_string_seams) { + // string NOT long enough to be worth aligning, skip + continue; + } - 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.); + // 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; }); - if (score < min_score) { - min_score = score; - out = sp.m_pos; + // repeat the alignment for the current seam, since it could be skipped due to alternative path being aligned. + global_index--; + + // gather all positions of seams and their weights + observations.resize(seam_string.size()); + observation_points.resize(seam_string.size()); + weights.resize(seam_string.size()); + + auto angle_3d = [](const Vec3f &a, const Vec3f &b) { return std::abs(acosf(a.normalized().dot(b.normalized()))); }; + + auto angle_weight = [](float angle) { return 1.0f / (0.1f + compute_angle_penalty(angle)); }; + + // gather points positions and weights + float total_length = 0.0f; + Vec3f last_point_pos = layers[seam_string[0].first].points[seam_string[0].second].position; + for (size_t index = 0; index < seam_string.size(); ++index) { + const SeamCandidate ¤t = layers[seam_string[index].first].points[seam_string[index].second]; + float layer_angle = 0.0f; + if (index > 0 && index < seam_string.size() - 1) { + layer_angle = angle_3d(current.position - layers[seam_string[index - 1].first].points[seam_string[index - 1].second].position, + layers[seam_string[index + 1].first].points[seam_string[index + 1].second].position - current.position); + } + observations[index] = current.position.head<2>(); + observation_points[index] = current.position.z(); + weights[index] = angle_weight(current.local_ccw_angle); + float sign = layer_angle > 2.0 * std::abs(current.local_ccw_angle) ? -1.0f : 1.0f; + if (current.type == EnforcedBlockedSeamPoint::Enforced) { + sign = 1.0f; + weights[index] += 3.0f; + } + total_length += sign * (last_point_pos - current.position).norm(); + last_point_pos = current.position; + } + + // Curve Fitting + size_t number_of_segments = std::max(size_t(1), size_t(std::max(0.0f, total_length) / SeamPlacer::seam_align_mm_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]; + float t = std::min(1.0f, std::abs(layers[pair.first].points[pair.second].local_ccw_angle) / SeamPlacer::sharp_angle_snapping_threshold); + if (layers[pair.first].points[pair.second].type == EnforcedBlockedSeamPoint::Enforced) { t = std::max(0.7f, t); } + + 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.0f - 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 } } - return out; +#ifdef DEBUG_FILES + fclose(clusters); + fclose(aligns); +#endif } - - -void SeamHistory::add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb) +void SeamPlacer::init(const Print &print, std::function throw_if_canceled_func) { - m_data_this_layer[po].push_back({pos, island_bb});; + using namespace SeamPlacerImpl; + m_seam_per_object.clear(); + + for (const PrintObject *po : print.objects()) { + throw_if_canceled_func(); + SeamPosition configured_seam_preference = po->config().seam_position.value; + SeamComparator comparator{configured_seam_preference}; + + { + GlobalModelInfo global_model_info{}; + gather_enforcers_blockers(global_model_info, po); + throw_if_canceled_func(); + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { compute_global_occlusion(global_model_info, po, throw_if_canceled_func); } + throw_if_canceled_func(); + 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"; + throw_if_canceled_func(); + 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"; + } + } // destruction of global_model_info (large structure, no longer needed) + throw_if_canceled_func(); + 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"; + throw_if_canceled_func(); + 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) + 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"; + } + throw_if_canceled_func(); + if (configured_seam_preference == spAligned || configured_seam_preference == spRear) { + 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(m_seam_per_object[po].layers, po->bounding_box(), comparator); +#endif + } } - - -void SeamHistory::clear() +void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const { - m_layer_id = 0; - m_data_last_layer.clear(); - m_data_this_layer.clear(); + 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))); + } + + 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 - 1 : seam_index - 1; + size_t index_of_next = seam_index == perimeter_point.perimeter.end_index - 1 ? perimeter_point.perimeter.start_index : seam_index + 1; + + 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; + + ExtrusionLoop::ClosestPathPoint 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.foot_pt)).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()); + } + + // Because the G-code export has 1um resolution, don't generate segments shorter than 1.5 microns, + // thus empty path segments will not be produced by G-code export. + //if (!loop.split_at_vertex(seam_point, scaled(0.0015))) { + 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 57c3532c3b..de99e12d34 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; } - - -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(); - -private: - struct SeamPoint { - Point m_pos; - BoundingBox m_island_bb; - }; - - std::map> m_data_last_layer; - std::map> m_data_this_layer; - size_t m_layer_id; -}; - - - -class SeamPlacer { -public: - void init(const Print& print); - - // 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. - - - 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); - - void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter, - const EdgeGrid::Grid* lower_layer_edge_grid); - - - using TreeType = AABBTreeIndirect::Tree<2, coord_t>; - using AlignedBoxType = Eigen::AlignedBox; - -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); - - 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()); - } -}; - +namespace EdgeGrid { +class Grid; } +namespace SeamPlacerImpl { + +// ************ FOR BACKPORT COMPATIBILITY ONLY *************** +// 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)); +} +// *************************** + +struct GlobalModelInfo; +struct SeamComparator; + +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 = Vec3f::Zero(); +}; + +// 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 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: + // Number of samples generated on the mesh. There are sqr_rays_per_sample_point*sqr_rays_per_sample_point rays casted from each samples + static constexpr size_t raycasting_visibility_samples_count = 30000; + // square of number of rays per sample point + static constexpr size_t sqr_rays_per_sample_point = 5; + + // arm length used during angles computation + static constexpr float polygon_local_angles_arm_distance = 0.3f; + // snapping angle - angles larger than this value will be snapped to during seam painting + static constexpr float sharp_angle_snapping_threshold = 55.0f * float(PI) / 180.0f; + // overhang angle for seam placement that still yields good results, in degrees, measured from vertical direction + static constexpr float overhang_angle_threshold = 45.0f * float(PI) / 180.0f; + + // determines angle importance compared to visibility ( neutral value is 1.0f. ) + static constexpr float angle_importance_aligned = 0.6f; + static constexpr float angle_importance_nearest = 1.0f; // use much higher angle importance for nearest mode, to combat the visibility info noise + + // 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; + + // 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.3f; + // seam_align_tolerable_dist_factor - how far to search for seam from current position, final dist is seam_align_tolerable_dist_factor * flow_width + static constexpr float seam_align_tolerable_dist_factor = 4.0f; + // minimum number of seams needed in cluster to make alignment happen + static constexpr size_t seam_align_minimum_string_seams = 6; + // millimeters covered by spline; determines number of splines for the given string + static constexpr size_t seam_align_mm_per_segment = 4.0f; + + // The following data structures hold all perimeter points for all PrintObject. + std::unordered_map m_seam_per_object; + + void init(const Print &print, std::function throw_if_canceled_func); + + void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const; + +private: + 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); + std::vector> find_seam_string(const PrintObject *po, std::pair start_seam, const SeamPlacerImpl::SeamComparator &comparator) const; + std::optional> find_next_seam_in_layer(const std::vector &layers, + const Vec3f & projected_position, + const size_t layer_idx, + const float max_distance, + const SeamPlacerImpl::SeamComparator & comparator) const; +}; + +} // namespace Slic3r + #endif // libslic3r_SeamPlacer_hpp_ diff --git a/src/libslic3r/Geometry/Bicubic.hpp b/src/libslic3r/Geometry/Bicubic.hpp new file mode 100644 index 0000000000..98c6f8bb29 --- /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/Curves.hpp b/src/libslic3r/Geometry/Curves.hpp new file mode 100644 index 0000000000..f4a5a0067e --- /dev/null +++ b/src/libslic3r/Geometry/Curves.hpp @@ -0,0 +1,218 @@ +#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; + + Vec get_fitted_value(const NumberType& value) const { + Vec 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 - int(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_linear_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); +} + +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 12e462569e..37c10827b1 100644 --- a/src/libslic3r/KDTreeIndirect.hpp +++ b/src/libslic3r/KDTreeIndirect.hpp @@ -11,224 +11,362 @@ namespace Slic3r { +enum class VisitorReturnMask : unsigned int { + CONTINUE_LEFT = 1, + CONTINUE_RIGHT = 2, + STOP = 4, +}; + // KD tree for N-dimensional closest point search. template class KDTreeIndirect { public: - static constexpr size_t NumDimensions = ANumDimensions; - using CoordinateFn = ACoordinateFn; - using CoordType = ACoordType; + static constexpr size_t NumDimensions = ANumDimensions; + using CoordinateFn = ACoordinateFn; + using CoordType = ACoordType; // Following could be static constexpr size_t, but that would not link in C++11 enum : size_t { npos = size_t(-1) }; - KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {} - KDTreeIndirect(CoordinateFn coordinate, std::vector indices) : coordinate(coordinate) { this->build(std::move(indices)); } - KDTreeIndirect(CoordinateFn coordinate, std::vector &&indices) : coordinate(coordinate) { this->build(std::move(indices)); } - KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); } - KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {} - KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; } - void clear() { m_nodes.clear(); } + KDTreeIndirect(CoordinateFn coordinate) : coordinate(coordinate) {} + KDTreeIndirect(CoordinateFn coordinate, std::vector indices) : coordinate(coordinate) { this->build(indices); } + KDTreeIndirect(CoordinateFn coordinate, size_t num_indices) : coordinate(coordinate) { this->build(num_indices); } + KDTreeIndirect(KDTreeIndirect &&rhs) : m_nodes(std::move(rhs.m_nodes)), coordinate(std::move(rhs.coordinate)) {} + KDTreeIndirect& operator=(KDTreeIndirect &&rhs) { m_nodes = std::move(rhs.m_nodes); coordinate = std::move(rhs.coordinate); return *this; } + void clear() { m_nodes.clear(); } - void build(size_t num_indices) - { - std::vector indices; - indices.reserve(num_indices); - for (size_t i = 0; i < num_indices; ++ i) - indices.emplace_back(i); - this->build(std::move(indices)); - } + void build(size_t num_indices) + { + std::vector indices; + indices.reserve(num_indices); + for (size_t i = 0; i < num_indices; ++ i) + indices.emplace_back(i); + this->build(indices); + } - void build(std::vector &&indices) - { - if (indices.empty()) - clear(); - else { - // Allocate enough memory for a full binary tree. - m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos); - build_recursive(indices, 0, 0, 0, indices.size() - 1); - } - indices.clear(); - } + void build(std::vector &indices) + { + if (indices.empty()) + clear(); + else { + // Allocate enough memory for a full binary tree. + m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos); + build_recursive(indices, 0, 0, 0, indices.size() - 1); + } + indices.clear(); + } - enum class VisitorReturnMask : unsigned int - { - CONTINUE_LEFT = 1, - CONTINUE_RIGHT = 2, - STOP = 4, - }; - 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); - return (dist * dist < search_radius + CoordType(EPSILON)) ? - // The plane intersects a hypersphere centered at point_coord of search_radius. - ((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) : - // The plane does not intersect the hypersphere. - (dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT); - } + 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); + return (dist * dist < search_radius + CoordType(EPSILON)) ? + // The plane intersects a hypersphere centered at point_coord of search_radius. + ((unsigned int)(VisitorReturnMask::CONTINUE_LEFT) | (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT)) : + // The plane does not intersect the hypersphere. + (dist > CoordType(0)) ? (unsigned int)(VisitorReturnMask::CONTINUE_RIGHT) : (unsigned int)(VisitorReturnMask::CONTINUE_LEFT); + } - // Visitor is supposed to return a bit mask of VisitorReturnMask. - template - void visit(Visitor &visitor) const - { + // Visitor is supposed to return a bit mask of VisitorReturnMask. + template + void visit(Visitor &visitor) const + { visit_recursive(0, 0, visitor); - } + } - CoordinateFn coordinate; + CoordinateFn coordinate; private: - // Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension. - void build_recursive(std::vector &input, size_t node, const size_t dimension, const size_t left, const size_t right) - { - if (left > right) - return; + // Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension. + void build_recursive(std::vector &input, size_t node, const size_t dimension, const size_t left, const size_t right) + { + if (left > right) + return; - assert(node < m_nodes.size()); + assert(node < m_nodes.size()); - if (left == right) { - // Insert a node into the balanced tree. - m_nodes[node] = input[left]; - return; - } + if (left == right) { + // Insert a node into the balanced tree. + m_nodes[node] = input[left]; + return; + } - // Partition the input to left / right pieces of the same length to produce a balanced tree. - size_t center = (left + right) / 2; - partition_input(input, dimension, left, right, center); - // Insert a node into the tree. - m_nodes[node] = input[center]; - // Build up the left / right subtrees. - size_t next_dimension = dimension; - if (++ next_dimension == NumDimensions) - next_dimension = 0; - if (center > left) - build_recursive(input, node * 2 + 1, next_dimension, left, center - 1); - build_recursive(input, node * 2 + 2, next_dimension, center + 1, right); - } + // Partition the input to left / right pieces of the same length to produce a balanced tree. + size_t center = (left + right) / 2; + partition_input(input, dimension, left, right, center); + // Insert a node into the tree. + m_nodes[node] = input[center]; + // Build up the left / right subtrees. + size_t next_dimension = dimension; + if (++ next_dimension == NumDimensions) + next_dimension = 0; + if (center > left) + build_recursive(input, node * 2 + 1, next_dimension, left, center - 1); + build_recursive(input, node * 2 + 2, next_dimension, center + 1, right); + } - // 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", - void partition_input(std::vector &input, const size_t dimension, size_t left, size_t right, const size_t k) const - { - while (left < right) { - size_t center = (left + right) / 2; - CoordType pivot; - { - // Bubble sort the input[left], input[center], input[right], so that a median of the three values - // will end up in input[center]. - CoordType left_value = this->coordinate(input[left], dimension); - CoordType center_value = this->coordinate(input[center], dimension); - CoordType right_value = this->coordinate(input[right], dimension); - if (left_value > center_value) { - std::swap(input[left], input[center]); - std::swap(left_value, center_value); - } - if (left_value > right_value) { - std::swap(input[left], input[right]); - right_value = left_value; - } - if (center_value > right_value) { - std::swap(input[center], input[right]); - center_value = right_value; - } - pivot = center_value; - } - if (right <= left + 2) - // The interval is already sorted. - break; - size_t i = left; - size_t j = right - 1; - std::swap(input[center], input[j]); - // Partition the set based on the pivot. - for (;;) { - // Skip left points that are already at correct positions. - // Search will certainly stop at position (right - 1), which stores the pivot. - while (this->coordinate(input[++ i], dimension) < pivot) ; - // Skip right points that are already at correct positions. - while (this->coordinate(input[-- j], dimension) > pivot && i < j) ; - if (i >= j) - break; - std::swap(input[i], input[j]); - } - // Restore pivot to the center of the sequence. - std::swap(input[i], input[right - 1]); - // Which side the kth element is in? - if (k < i) - right = i - 1; - else if (k == i) - // Sequence is partitioned, kth element is at its place. - break; - else - left = i + 1; - } - } + // 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", + void partition_input(std::vector &input, const size_t dimension, size_t left, size_t right, const size_t k) const + { + while (left < right) { + size_t center = (left + right) / 2; + CoordType pivot; + { + // Bubble sort the input[left], input[center], input[right], so that a median of the three values + // will end up in input[center]. + CoordType left_value = this->coordinate(input[left], dimension); + CoordType center_value = this->coordinate(input[center], dimension); + CoordType right_value = this->coordinate(input[right], dimension); + if (left_value > center_value) { + std::swap(input[left], input[center]); + std::swap(left_value, center_value); + } + if (left_value > right_value) { + std::swap(input[left], input[right]); + right_value = left_value; + } + if (center_value > right_value) { + std::swap(input[center], input[right]); + center_value = right_value; + } + pivot = center_value; + } + if (right <= left + 2) + // The interval is already sorted. + break; + size_t i = left; + size_t j = right - 1; + std::swap(input[center], input[j]); + // Partition the set based on the pivot. + for (;;) { + // Skip left points that are already at correct positions. + // Search will certainly stop at position (right - 1), which stores the pivot. + while (this->coordinate(input[++ i], dimension) < pivot) ; + // Skip right points that are already at correct positions. + while (this->coordinate(input[-- j], dimension) > pivot && i < j) ; + if (i >= j) + break; + std::swap(input[i], input[j]); + } + // Restore pivot to the center of the sequence. + std::swap(input[i], input[right - 1]); + // Which side the kth element is in? + if (k < i) + right = i - 1; + else if (k == i) + // Sequence is partitioned, kth element is at its place. + break; + else + left = i + 1; + } + } - template - void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const - { - assert(! m_nodes.empty()); - if (node >= m_nodes.size() || m_nodes[node] == npos) - return; + template + void visit_recursive(size_t node, size_t dimension, Visitor &visitor) const + { + assert(! m_nodes.empty()); + if (node >= m_nodes.size() || m_nodes[node] == npos) + return; - // Left / right child node index. - size_t left = node * 2 + 1; - size_t right = left + 1; - unsigned int mask = visitor(m_nodes[node], dimension); - if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) { - size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension; - if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT) - visit_recursive(left, next_dimension, visitor); - if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT) - visit_recursive(right, next_dimension, visitor); - } - } + // Left / right child node index. + size_t left = node * 2 + 1; + size_t right = left + 1; + unsigned int mask = visitor(m_nodes[node], dimension); + if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) { + size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension; + if (mask & (unsigned int)VisitorReturnMask::CONTINUE_LEFT) + visit_recursive(left, next_dimension, visitor); + if (mask & (unsigned int)VisitorReturnMask::CONTINUE_RIGHT) + visit_recursive(right, next_dimension, visitor); + } + } - std::vector m_nodes; + std::vector m_nodes; }; // Find a closest point using Euclidian metrics. // Returns npos if not found. -template -size_t find_closest_point(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) +template +std::array find_closest_points( + const KDTreeIndirect &kdtree, + const PointType &point, + FilterFn filter) { - using CoordType = typename KDTreeIndirectType::CoordType; + using Tree = KDTreeIndirect; - struct Visitor { - const KDTreeIndirectType &kdtree; - const PointType &point; - const FilterFn filter; - size_t min_idx = KDTreeIndirectType::npos; - CoordType min_dist = std::numeric_limits::max(); + struct Visitor + { + const Tree &kdtree; + const PointType &point; + const FilterFn filter; - Visitor(const KDTreeIndirectType &kdtree, const PointType &point, FilterFn filter) : kdtree(kdtree), point(point), 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 = point[i] - kdtree.coordinate(idx, i); - dist += d * d; - } - if (dist < min_dist) { - min_dist = dist; - min_idx = idx; - } - } - return kdtree.descent_mask(point[dimension], min_dist, idx, dimension); - } - } visitor(kdtree, point, filter); + std::array, K> results; - kdtree.visit(visitor); - return visitor.min_idx; + Visitor(const Tree &kdtree, const PointType &point, FilterFn filter) + : kdtree(kdtree), point(point), filter(filter) + { + results.fill(std::make_pair(Tree::npos, + std::numeric_limits::max())); + } + unsigned int operator()(size_t idx, size_t dimension) + { + if (this->filter(idx)) { + auto dist = CoordT(0); + for (size_t i = 0; i < D; ++i) { + CoordT d = point[i] - kdtree.coordinate(idx, i); + dist += d * d; + } + + auto res = std::make_pair(idx, dist); + auto it = std::lower_bound(results.begin(), results.end(), + res, [](auto &r1, auto &r2) { + return r1.second < r2.second; + }); + + if (it != results.end()) { + std::rotate(it, std::prev(results.end()), results.end()); + *it = res; + } + } + return kdtree.descent_mask(point[dimension], + results.front().second, idx, + dimension); + } + } visitor(kdtree, point, filter); + + kdtree.visit(visitor); + std::array ret; + for (size_t i = 0; i < K; i++) ret[i] = visitor.results[i].first; + + return ret; +} + +template +std::array find_closest_points( + const KDTreeIndirect &kdtree, const PointType &point) +{ + return find_closest_points(kdtree, point, [](size_t) { return true; }); +} + +template +size_t find_closest_point(const KDTreeIndirect &kdtree, + const PointType &point, + FilterFn filter) +{ + return find_closest_points<1>(kdtree, point, filter)[0]; } template size_t find_closest_point(const KDTreeIndirectType& kdtree, const PointType& point) { - return find_closest_point(kdtree, point, [](size_t) { return true; }); + 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; + }); +} + +// Find nearby points (spherical neighbourhood) using Euclidian metrics. +template +std::vector find_nearby_points(const KDTreeIndirectType &kdtree, + const PointType &bb_min, + const PointType &bb_max, + FilterFn filter) +{ + struct Visitor { + const KDTreeIndirectType &kdtree; + const PointType &bb_min, &bb_max; + const FilterFn filter; + std::vector result; + + Visitor(const KDTreeIndirectType &kdtree, const PointType& bbmin, const PointType& bbmax, + FilterFn filter) : + kdtree(kdtree), bb_min{bbmin}, bb_max{bbmax}, filter(filter) { + } + unsigned int operator()(size_t idx, size_t dimension) { + unsigned int ret = + static_cast(VisitorReturnMask::CONTINUE_LEFT) | + static_cast(VisitorReturnMask::CONTINUE_RIGHT); + + if (this->filter(idx)) { + PointType p; + bool contains = true; + for (size_t i = 0; i < KDTreeIndirectType::NumDimensions; ++i) { + p(i) = kdtree.coordinate(idx, i); + contains = contains && bb_min(i) <= p(i) && p(i) <= bb_max(i); + } + + if (p(dimension) < bb_min(dimension)) + ret = static_cast(VisitorReturnMask::CONTINUE_RIGHT); + if (p(dimension) > bb_max(dimension)) + ret = static_cast(VisitorReturnMask::CONTINUE_LEFT); + + if (contains) + result.emplace_back(idx); + } + + return ret; + } + } visitor(kdtree, bb_min, bb_max, filter); + + kdtree.visit(visitor); + return visitor.result; } } // namespace Slic3r diff --git a/src/libslic3r/MultiPoint.cpp b/src/libslic3r/MultiPoint.cpp index 227f0eea3d..0076300ea6 100644 --- a/src/libslic3r/MultiPoint.cpp +++ b/src/libslic3r/MultiPoint.cpp @@ -63,6 +63,23 @@ int MultiPoint::find_point(const Point &point) const return -1; // not found } +int MultiPoint::find_point(const Point &point, double scaled_epsilon) const +{ + if (scaled_epsilon == 0) return this->find_point(point); + + auto dist2_min = std::numeric_limits::max(); + auto eps2 = scaled_epsilon * scaled_epsilon; + int idx_min = -1; + for (const Point &pt : this->points) { + double d2 = (pt - point).cast().squaredNorm(); + if (d2 < dist2_min) { + idx_min = int(&pt - &this->points.front()); + dist2_min = d2; + } + } + return dist2_min < eps2 ? idx_min : -1; +} + bool MultiPoint::has_boundary_point(const Point &point) const { double dist = (point.projection_onto(*this) - point).cast().norm(); diff --git a/src/libslic3r/MultiPoint.hpp b/src/libslic3r/MultiPoint.hpp index 86555c33a6..4c1f9049ed 100644 --- a/src/libslic3r/MultiPoint.hpp +++ b/src/libslic3r/MultiPoint.hpp @@ -43,7 +43,12 @@ public: double length() const; bool is_valid() const { return this->points.size() >= 2; } + // Return index of a polygon point exactly equal to point. + // Return -1 if no such point exists. int find_point(const Point &point) const; + // Return index of the closest point to point closer than scaled_epsilon. + // Return -1 if no such point exists. + int find_point(const Point &point, const double scaled_epsilon) const; bool has_boundary_point(const Point &point) const; int closest_point_index(const Point &point) const { int idx = -1; diff --git a/src/libslic3r/NormalUtils.cpp b/src/libslic3r/NormalUtils.cpp new file mode 100644 index 0000000000..dc94515656 --- /dev/null +++ b/src/libslic3r/NormalUtils.cpp @@ -0,0 +1,142 @@ +#include "NormalUtils.hpp" + +using namespace Slic3r; + +Vec3f NormalUtils::create_triangle_normal( + const stl_triangle_vertex_indices &indices, + const std::vector & vertices) +{ + const stl_vertex &v0 = vertices[indices[0]]; + const stl_vertex &v1 = vertices[indices[1]]; + const stl_vertex &v2 = vertices[indices[2]]; + Vec3f direction = (v1 - v0).cross(v2 - v0); + direction.normalize(); + return direction; +} + +std::vector NormalUtils::create_triangle_normals( + const indexed_triangle_set &its) +{ + std::vector normals; + normals.reserve(its.indices.size()); + for (const Vec3crd &index : its.indices) { + normals.push_back(create_triangle_normal(index, its.vertices)); + } + return normals; +} + +NormalUtils::Normals NormalUtils::create_normals_average_neighbor( + const indexed_triangle_set &its) +{ + size_t count_vertices = its.vertices.size(); + std::vector normals(count_vertices, Vec3f(.0f, .0f, .0f)); + std::vector count(count_vertices, 0); + for (const Vec3crd &indice : its.indices) { + Vec3f normal = create_triangle_normal(indice, its.vertices); + for (int i = 0; i < 3; ++i) { + normals[indice[i]] += normal; + ++count[indice[i]]; + } + } + // normalize to size 1 + for (auto &normal : normals) { + size_t index = &normal - &normals.front(); + normal /= static_cast(count[index]); + } + return normals; +} + +// calc triangle angle of vertex defined by index to triangle indices +float NormalUtils::indice_angle(int i, + const Vec3crd & indice, + const std::vector &vertices) +{ + int i1 = (i == 0) ? 2 : (i - 1); + int i2 = (i == 2) ? 0 : (i + 1); + + Vec3f v1 = vertices[i1] - vertices[i]; + Vec3f v2 = vertices[i2] - vertices[i]; + + v1.normalize(); + v2.normalize(); + + float w = v1.dot(v2); + if (w > 1.f) + w = 1.f; + else if (w < -1.f) + w = -1.f; + return acos(w); +} + +NormalUtils::Normals NormalUtils::create_normals_angle_weighted( + const indexed_triangle_set &its) +{ + size_t count_vertices = its.vertices.size(); + std::vector normals(count_vertices, Vec3f(.0f, .0f, .0f)); + std::vector count(count_vertices, 0.f); + for (const Vec3crd &indice : its.indices) { + Vec3f normal = create_triangle_normal(indice, its.vertices); + Vec3f angles(indice_angle(0, indice, its.vertices), + indice_angle(1, indice, its.vertices), 0.f); + angles[2] = (M_PI - angles[0] - angles[1]); + for (int i = 0; i < 3; ++i) { + const float &weight = angles[i]; + normals[indice[i]] += normal * weight; + count[indice[i]] += weight; + } + } + // normalize to size 1 + for (auto &normal : normals) { + size_t index = &normal - &normals.front(); + normal /= count[index]; + } + return normals; +} + +NormalUtils::Normals NormalUtils::create_normals_nelson_weighted( + const indexed_triangle_set &its) +{ + size_t count_vertices = its.vertices.size(); + std::vector normals(count_vertices, Vec3f(.0f, .0f, .0f)); + std::vector count(count_vertices, 0.f); + const std::vector &vertices = its.vertices; + for (const Vec3crd &indice : its.indices) { + Vec3f normal = create_triangle_normal(indice, vertices); + + const stl_vertex &v0 = vertices[indice[0]]; + const stl_vertex &v1 = vertices[indice[1]]; + const stl_vertex &v2 = vertices[indice[2]]; + + float e0 = (v0 - v1).norm(); + float e1 = (v1 - v2).norm(); + float e2 = (v2 - v0).norm(); + + Vec3f coefs(e0 * e2, e0 * e1, e1 * e2); + for (int i = 0; i < 3; ++i) { + const float &weight = coefs[i]; + normals[indice[i]] += normal * weight; + count[indice[i]] += weight; + } + } + // normalize to size 1 + for (auto &normal : normals) { + size_t index = &normal - &normals.front(); + normal /= count[index]; + } + return normals; +} + +// calculate normals by averaging normals of neghbor triangles +std::vector NormalUtils::create_normals( + const indexed_triangle_set &its, VertexNormalType type) +{ + switch (type) { + case VertexNormalType::AverageNeighbor: + return create_normals_average_neighbor(its); + case VertexNormalType::AngleWeighted: + return create_normals_angle_weighted(its); + case VertexNormalType::NelsonMaxWeighted: + default: + return create_normals_nelson_weighted(its); + } +} diff --git a/src/libslic3r/NormalUtils.hpp b/src/libslic3r/NormalUtils.hpp new file mode 100644 index 0000000000..60ec57f724 --- /dev/null +++ b/src/libslic3r/NormalUtils.hpp @@ -0,0 +1,69 @@ +#ifndef slic3r_NormalUtils_hpp_ +#define slic3r_NormalUtils_hpp_ + +#include "Point.hpp" +#include "Model.hpp" + +namespace Slic3r { + +/// +/// Collection of static function +/// to create normals +/// +class NormalUtils +{ +public: + using Normal = Vec3f; + using Normals = std::vector; + NormalUtils() = delete; // only static functions + + enum class VertexNormalType { + AverageNeighbor, + AngleWeighted, + NelsonMaxWeighted + }; + + /// + /// Create normal for triangle defined by indices from vertices + /// + /// index into vertices + /// vector of vertices + /// normal to triangle(normalized to size 1) + static Normal create_triangle_normal( + const stl_triangle_vertex_indices &indices, + const std::vector & vertices); + + /// + /// Create normals for each vertices + /// + /// indices and vertices + /// Vector of normals + static Normals create_triangle_normals(const indexed_triangle_set &its); + + /// + /// Create normals for each vertex by averaging neighbor triangles normal + /// + /// Triangle indices and vertices + /// Type of calculation normals + /// Normal for each vertex + static Normals create_normals( + const indexed_triangle_set &its, + VertexNormalType type = VertexNormalType::NelsonMaxWeighted); + static Normals create_normals_average_neighbor(const indexed_triangle_set &its); + static Normals create_normals_angle_weighted(const indexed_triangle_set &its); + static Normals create_normals_nelson_weighted(const indexed_triangle_set &its); + + /// + /// Calculate angle of trinagle side. + /// + /// index to indices, define angle point + /// address to vertices + /// vertices data + /// Angle [in radian] + static float indice_angle(int i, + const Vec3crd & indice, + const std::vector &vertices); +}; + +} // namespace Slic3r +#endif // slic3r_NormalUtils_hpp_ diff --git a/src/libslic3r/Polyline.cpp b/src/libslic3r/Polyline.cpp index ca11357bb4..78b064e24f 100644 --- a/src/libslic3r/Polyline.cpp +++ b/src/libslic3r/Polyline.cpp @@ -517,6 +517,28 @@ bool remove_degenerate(Polylines &polylines) return modified; } +std::pair foot_pt(const Points &polyline, const Point &pt) +{ + if (polyline.size() < 2) return std::make_pair(-1, Point(0, 0)); + + auto d2_min = std::numeric_limits::max(); + Point foot_pt_min; + Point prev = polyline.front(); + auto it = polyline.begin(); + auto it_proj = polyline.begin(); + for (++it; it != polyline.end(); ++it) { + Point foot_pt = pt.projection_onto(Line(prev, *it)); + double d2 = (foot_pt - pt).cast().squaredNorm(); + if (d2 < d2_min) { + d2_min = d2; + foot_pt_min = foot_pt; + it_proj = it; + } + prev = *it; + } + return std::make_pair(int(it_proj - polyline.begin()) - 1, foot_pt_min); +} + ThickLines ThickPolyline::thicklines() const { ThickLines lines; diff --git a/src/libslic3r/Polyline.hpp b/src/libslic3r/Polyline.hpp index d97801ac63..45361f68e2 100644 --- a/src/libslic3r/Polyline.hpp +++ b/src/libslic3r/Polyline.hpp @@ -222,6 +222,9 @@ const Point& leftmost_point(const Polylines &polylines); bool remove_degenerate(Polylines &polylines); +// Returns index of a segment of a polyline and foot point of pt on polyline. +std::pair foot_pt(const Points &polyline, const Point &pt); + class ThickPolyline : public Polyline { public: ThickPolyline() : endpoints(std::make_pair(false, false)) {} diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp new file mode 100644 index 0000000000..b36278c37f --- /dev/null +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -0,0 +1,183 @@ +#include "ShortEdgeCollapse.hpp" +#include "libslic3r/NormalUtils.hpp" + +#include +#include +#include +#include + +namespace Slic3r { + +void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) { + // whenever vertex is removed, its mapping is update to the index of vertex with wich it merged + std::vector vertices_index_mapping(mesh.vertices.size()); + for (size_t idx = 0; idx < vertices_index_mapping.size(); ++idx) { + vertices_index_mapping[idx] = idx; + } + // Algorithm uses get_final_index query to get the actual vertex index. The query also updates all mappings on the way, essentially flattening the mapping + std::vector flatten_queue; + auto get_final_index = [&vertices_index_mapping, &flatten_queue](const size_t &orig_index) { + flatten_queue.clear(); + size_t idx = orig_index; + while (vertices_index_mapping[idx] != idx) { + flatten_queue.push_back(idx); + idx = vertices_index_mapping[idx]; + } + for (size_t i : flatten_queue) { + vertices_index_mapping[i] = idx; + } + return idx; + + }; + + // if face is removed, mark it here + std::vector face_removal_flags(mesh.indices.size(), false); + + std::vector triangles_neighbors = its_face_neighbors_par(mesh); + + // now compute vertices dot product - this is used during edge collapse, + // to determine which vertex to remove and which to keep; We try to keep the one with larger angle, because it defines the shape "more". + // The min vertex dot product is lowest dot product of its normal with the normals of faces around it. + // the lower the dot product, the more we want to keep the vertex + // NOTE: This score is not updated, even though the decimation does change the mesh. It saves computation time, and there are no strong reasons to update. + std::vector min_vertex_dot_product(mesh.vertices.size(), 1); + { + std::vector face_normals = its_face_normals(mesh); + std::vector vertex_normals = NormalUtils::create_normals(mesh); + + for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) { + Vec3i t = mesh.indices[face_idx]; + Vec3f n = face_normals[face_idx]; + min_vertex_dot_product[t[0]] = std::min(min_vertex_dot_product[t[0]], n.dot(vertex_normals[t[0]])); + min_vertex_dot_product[t[1]] = std::min(min_vertex_dot_product[t[1]], n.dot(vertex_normals[t[1]])); + min_vertex_dot_product[t[2]] = std::min(min_vertex_dot_product[t[2]], n.dot(vertex_normals[t[2]])); + } + } + + // lambda to remove face. It flags the face as removed, and updates neighbourhood info + auto remove_face = [&triangles_neighbors, &face_removal_flags](int face_idx, int other_face_idx) { + if (face_idx < 0) { + return; + } + face_removal_flags[face_idx] = true; + Vec3i neighbors = triangles_neighbors[face_idx]; + int n_a = neighbors[0] != other_face_idx ? neighbors[0] : neighbors[1]; + int n_b = neighbors[2] != other_face_idx ? neighbors[2] : neighbors[1]; + if (n_a > 0) + for (int &n : triangles_neighbors[n_a]) { + if (n == face_idx) { + n = n_b; + break; + } + } + if (n_b > 0) + for (int &n : triangles_neighbors[n_b]) { + if (n == face_idx) { + n = n_a; + break; + } + } + }; + + std::mt19937_64 generator { 27644437 };// default constant seed! so that results are deterministic + std::vector face_indices(mesh.indices.size()); + for (size_t idx = 0; idx < face_indices.size(); ++idx) { + face_indices[idx] = idx; + } + //tmp face indices used only for swapping + std::vector tmp_face_indices(mesh.indices.size()); + + float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all + float edge_len = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased + + while (face_indices.size() > target_triangle_count) { + // simpple func to increase the edge len - if decimation ratio is low, it increases the len up to twice, if decimation ratio is high, increments are low + edge_len = edge_len * (1.0f + 1.0 - decimation_ratio); + float max_edge_len_squared = edge_len * edge_len; + + //shuffle the faces and traverse in random order, this MASSIVELY improves the quality of the result + std::shuffle(face_indices.begin(), face_indices.end(), generator); + + for (const size_t &face_idx : face_indices) { + if (face_removal_flags[face_idx]) { + // if face already removed from previous collapses, skip (each collapse removes two triangles [at least] ) + continue; + } + + // look at each edge if it is good candidate for collapse + for (size_t edge_idx = 0; edge_idx < 3; ++edge_idx) { + size_t vertex_index_keep = get_final_index(mesh.indices[face_idx][edge_idx]); + size_t vertex_index_remove = get_final_index(mesh.indices[face_idx][(edge_idx + 1) % 3]); + //check distance, skip long edges + if ((mesh.vertices[vertex_index_keep] - mesh.vertices[vertex_index_remove]).squaredNorm() + > max_edge_len_squared) { + continue; + } + // swap indexes if vertex_index_keep has higher dot product (we want to keep low dot product vertices) + if (min_vertex_dot_product[vertex_index_remove] < min_vertex_dot_product[vertex_index_keep]) { + size_t tmp = vertex_index_keep; + vertex_index_keep = vertex_index_remove; + vertex_index_remove = tmp; + } + + //remove vertex + { + // map its index to the index of the kept vertex + vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep]; + } + + int neighbor_to_remove_face_idx = triangles_neighbors[face_idx][edge_idx]; + // remove faces + remove_face(face_idx, neighbor_to_remove_face_idx); + remove_face(neighbor_to_remove_face_idx, face_idx); + + // break. this triangle is done + break; + } + } + + // filter face_indices, remove those that have been collapsed + size_t prev_size = face_indices.size(); + tmp_face_indices.clear(); + for (size_t face_idx : face_indices) { + if (!face_removal_flags[face_idx]){ + tmp_face_indices.push_back(face_idx); + } + } + face_indices.swap(tmp_face_indices); + + decimation_ratio = float(prev_size - face_indices.size()) / float(prev_size); + //std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; + } + + //Extract the result mesh + std::unordered_map final_vertices_mapping; + std::vector final_vertices; + std::vector final_indices; + final_indices.reserve(face_indices.size()); + for (size_t idx : face_indices) { + Vec3i final_face; + for (size_t i = 0; i < 3; ++i) { + final_face[i] = get_final_index(mesh.indices[idx][i]); + } + if (final_face[0] == final_face[1] || final_face[1] == final_face[2] || final_face[2] == final_face[0]) { + continue; // discard degenerate triangles + } + + for (size_t i = 0; i < 3; ++i) { + if (final_vertices_mapping.find(final_face[i]) == final_vertices_mapping.end()) { + final_vertices_mapping[final_face[i]] = final_vertices.size(); + final_vertices.push_back(mesh.vertices[final_face[i]]); + } + final_face[i] = final_vertices_mapping[final_face[i]]; + } + + final_indices.push_back(final_face); + } + + mesh.vertices = final_vertices; + mesh.indices = final_indices; +} + +} //namespace Slic3r + diff --git a/src/libslic3r/ShortEdgeCollapse.hpp b/src/libslic3r/ShortEdgeCollapse.hpp new file mode 100644 index 0000000000..e6f1822c89 --- /dev/null +++ b/src/libslic3r/ShortEdgeCollapse.hpp @@ -0,0 +1,16 @@ +#ifndef SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_ +#define SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_ + +#include "libslic3r/TriangleMesh.hpp" + +namespace Slic3r{ + +// Decimates the model by collapsing short edges. It starts with very small edges and gradually increases the collapsible length, +// until the target triangle count is reached (the algorithm will certainly undershoot the target count, result will have less triangles than target count) +// The algorithm does not check for triangle flipping, disconnections, self intersections or any other degeneration that can appear during mesh processing. +void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count); + +} + + +#endif /* SRC_LIBSLIC3R_SHORTEDGECOLLAPSE_HPP_ */ diff --git a/src/libslic3r/TriangleSetSampling.cpp b/src/libslic3r/TriangleSetSampling.cpp new file mode 100644 index 0000000000..bb03ff6d75 --- /dev/null +++ b/src/libslic3r/TriangleSetSampling.cpp @@ -0,0 +1,71 @@ +#include "TriangleSetSampling.hpp" +#include +#include +#include +#include + +namespace Slic3r { + +TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set) { + std::vector triangles_area(triangle_set.indices.size()); + + tbb::parallel_for(tbb::blocked_range(0, triangle_set.indices.size()), + [&triangle_set, &triangles_area]( + tbb::blocked_range r) { + for (size_t t_idx = r.begin(); t_idx < r.end(); ++t_idx) { + const Vec3f &a = triangle_set.vertices[triangle_set.indices[t_idx].x()]; + const Vec3f &b = triangle_set.vertices[triangle_set.indices[t_idx].y()]; + const Vec3f &c = triangle_set.vertices[triangle_set.indices[t_idx].z()]; + double area = double(0.5 * (b - a).cross(c - a).norm()); + triangles_area[t_idx] = area; + } + }); + + std::map area_sum_to_triangle_idx; + float area_sum = 0; + for (size_t t_idx = 0; t_idx < triangles_area.size(); ++t_idx) { + area_sum += triangles_area[t_idx]; + area_sum_to_triangle_idx[area_sum] = t_idx; + } + + std::mt19937_64 mersenne_engine { 27644437 }; + // random numbers on interval [0, 1) + std::uniform_real_distribution fdistribution; + + auto get_random = [&fdistribution, &mersenne_engine]() { + return Vec3d { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) }; + }; + + std::vector random_samples(samples_count); + std::generate(random_samples.begin(), random_samples.end(), get_random); + + TriangleSetSamples result; + result.total_area = area_sum; + result.positions.resize(samples_count); + result.normals.resize(samples_count); + result.triangle_indices.resize(samples_count); + + tbb::parallel_for(tbb::blocked_range(0, samples_count), + [&triangle_set, &area_sum_to_triangle_idx, &area_sum, &random_samples, &result]( + tbb::blocked_range r) { + for (size_t s_idx = r.begin(); s_idx < r.end(); ++s_idx) { + double t_sample = random_samples[s_idx].x() * area_sum; + size_t t_idx = area_sum_to_triangle_idx.upper_bound(t_sample)->second; + + double sq_u = std::sqrt(random_samples[s_idx].y()); + double v = random_samples[s_idx].z(); + + Vec3f A = triangle_set.vertices[triangle_set.indices[t_idx].x()]; + Vec3f B = triangle_set.vertices[triangle_set.indices[t_idx].y()]; + Vec3f C = triangle_set.vertices[triangle_set.indices[t_idx].z()]; + + result.positions[s_idx] = A * (1 - sq_u) + B * (sq_u * (1 - v)) + C * (v * sq_u); + result.normals[s_idx] = ((B - A).cross(C - B)).normalized(); + result.triangle_indices[s_idx] = t_idx; + } + }); + + return result; +} + +} diff --git a/src/libslic3r/TriangleSetSampling.hpp b/src/libslic3r/TriangleSetSampling.hpp new file mode 100644 index 0000000000..28a661d76f --- /dev/null +++ b/src/libslic3r/TriangleSetSampling.hpp @@ -0,0 +1,20 @@ +#ifndef SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_ +#define SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_ + +#include +#include "libslic3r/Point.hpp" + +namespace Slic3r { + +struct TriangleSetSamples { + float total_area; + std::vector positions; + std::vector normals; + std::vector triangle_indices; +}; + +TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const indexed_triangle_set &triangle_set); + +} + +#endif /* SRC_LIBSLIC3R_TRIANGLESETSAMPLING_HPP_ */ diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 63b0422907..44d5de1305 100644 --- a/src/libslic3r/libslic3r.h +++ b/src/libslic3r/libslic3r.h @@ -24,6 +24,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" @@ -96,6 +103,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); } @@ -341,6 +358,12 @@ public: inline bool empty() const { return size() == 0; } }; +template> +constexpr T NaN = std::numeric_limits::quiet_NaN(); + +constexpr float NaNf = NaN; +constexpr double NaNd = NaN; + } // namespace Slic3r #endif