From d59d8eebdea068ae0ce1efe8a9c464a9c1237cb3 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Fri, 20 May 2022 14:47:05 +0200 Subject: [PATCH 01/11] Full occlusion guided seam placer backport --- src/libslic3r/AABBTreeIndirect.hpp | 160 +- src/libslic3r/AABBTreeLines.hpp | 112 ++ src/libslic3r/CMakeLists.txt | 4 + src/libslic3r/GCode.cpp | 62 +- src/libslic3r/GCode.hpp | 6 +- src/libslic3r/GCode/SeamPlacer.cpp | 2343 +++++++++++++++---------- src/libslic3r/GCode/SeamPlacer.hpp | 248 ++- src/libslic3r/Geometry/Bicubic.hpp | 291 +++ src/libslic3r/Geometry/Curves.hpp | 205 +++ src/libslic3r/KDTreeIndirect.hpp | 498 ++++-- src/libslic3r/QuadricEdgeCollapse.cpp | 14 +- src/libslic3r/SLA/bicubic.h | 186 -- src/libslic3r/Subdivide.cpp | 218 +++ src/libslic3r/Subdivide.hpp | 12 + src/libslic3r/libslic3r.h | 23 + 15 files changed, 2889 insertions(+), 1493 deletions(-) create mode 100644 src/libslic3r/AABBTreeLines.hpp create mode 100644 src/libslic3r/Geometry/Bicubic.hpp create mode 100644 src/libslic3r/Geometry/Curves.hpp delete mode 100644 src/libslic3r/SLA/bicubic.h create mode 100644 src/libslic3r/Subdivide.cpp create mode 100644 src/libslic3r/Subdivide.hpp diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp index 217166f8c..9b9c886ec 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 000000000..7b9595419 --- /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 16597e649..21fac7456 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -128,10 +128,12 @@ add_library(libslic3r STATIC 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 @@ -215,6 +217,8 @@ add_library(libslic3r STATIC SlicesToTriangleMesh.cpp SlicingAdaptive.cpp SlicingAdaptive.hpp + Subdivide.cpp + Subdivide.hpp SupportMaterial.cpp SupportMaterial.hpp Surface.cpp diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index d3310fe70..d91078164 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -3,7 +3,6 @@ #include "GCode.hpp" #include "Exception.hpp" #include "ExtrusionEntity.hpp" -#include "EdgeGrid.hpp" #include "Geometry/ConvexHull.hpp" #include "GCode/PrintExtents.hpp" #include "GCode/WipeTower.hpp" @@ -2395,7 +2394,6 @@ GCode::LayerResult GCode::process_layer( } // for objects // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders. - std::vector> lower_layer_edge_grids(layers.size()); for (unsigned int extruder_id : layer_tools.extruders) { gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ? @@ -2489,9 +2487,9 @@ GCode::LayerResult GCode::process_layer( //FIXME the following code prints regions in the order they are defined, the path is not optimized in any way. if (print.config().infill_first) { gcode += this->extrude_infill(print, by_region_specific, false); - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); } else { - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); gcode += this->extrude_infill(print,by_region_specific, false); } // ironing @@ -2626,51 +2624,22 @@ std::string GCode::change_layer(coordf_t print_z) return gcode; } - - -static std::unique_ptr calculate_layer_edge_grid(const Layer& layer) -{ - auto out = make_unique(); - - // Create the distance field for a layer below. - const coord_t distance_field_resolution = coord_t(scale_(1.) + 0.5); - out->create(layer.lslices, distance_field_resolution); - out->calculate_sdf(); -#if 0 - { - static int iRun = 0; - BoundingBox bbox = (*lower_layer_edge_grid)->bbox(); - bbox.min(0) -= scale_(5.f); - bbox.min(1) -= scale_(5.f); - bbox.max(0) += scale_(5.f); - bbox.max(1) += scale_(5.f); - EdgeGrid::save_png(*(*lower_layer_edge_grid), bbox, scale_(0.1f), debug_out_path("GCode_extrude_loop_edge_grid-%d.png", iRun++)); - } -#endif - return out; -} - - -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); - // extrude all loops ccw bool was_clockwise = loop.make_counter_clockwise(); // find the point of the loop that is closest to the current extruder position // or randomize if requested Point last_pos = this->last_pos(); - if (m_config.spiral_vase) { + if (! m_config.spiral_vase && description == "perimeter") { + assert(m_layer != nullptr); + m_seam_placer.place_seam(m_layer, loop, m_config.external_perimeters_first, this->last_pos()); + } else loop.split_at(last_pos, false); - } - else - m_seam_placer.place_seam(loop, this->last_pos(), m_config.external_perimeters_first, - EXTRUDER_CONFIG(nozzle_diameter), lower_layer_edge_grid ? lower_layer_edge_grid->get() : nullptr); // clip the path to avoid the extruder to get exactly on the first point of the loop; // if polyline was shorter than the clipping distance we'd get a null polyline, so @@ -2759,14 +2728,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 ""; @@ -2787,24 +2756,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 e37b47d81..08b212007 100644 --- a/src/libslic3r/GCode.hpp +++ b/src/libslic3r/GCode.hpp @@ -271,8 +271,8 @@ private: void set_extruders(const std::vector &extruder_ids); std::string preamble(); std::string change_layer(coordf_t print_z); - std::string extrude_entity(const ExtrusionEntity &entity, std::string description = "", double speed = -1., std::unique_ptr *lower_layer_edge_grid = nullptr); - std::string extrude_loop(ExtrusionLoop loop, std::string description, double speed = -1., std::unique_ptr *lower_layer_edge_grid = nullptr); + std::string extrude_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.); @@ -339,7 +339,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 bcd238a72..4f30fdc19 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -1,1011 +1,1550 @@ #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/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/QuadricEdgeCollapse.hpp" +#include "libslic3r/Subdivide.hpp" + +#include "libslic3r/Geometry/Curves.hpp" + +#include "libslic3r/Utils.hpp" + +//#define DEBUG_FILES + +#ifdef DEBUG_FILES +#include +#include +#endif namespace Slic3r { -// This penalty is added to all points inside custom blockers (subtracted from pts inside enforcers). -static constexpr float ENFORCER_BLOCKER_PENALTY = 100; +namespace SeamPlacerImpl { -// In case there are custom enforcers/blockers, the loop polygon shall always have -// sides smaller than this (so it isn't limited to original resolution). -static constexpr float MINIMAL_POLYGON_SIDE = scaled(0.2f); - -// When spAligned is active and there is a support enforcer, -// add this penalty to its center. -static constexpr float ENFORCER_CENTER_PENALTY = -10.f; - - - -// This function was introduced in 2016 to assign penalties to overhangs. -// LukasM thinks that it discriminated a bit too much, so especially external -// seams were than placed in funny places (non-overhangs were preferred too much). -// He implemented his own version (below) which applies fixed penalty for really big overlaps. -// static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) -// { -// // The extrudate is not fully supported by the lower layer. Fit a polynomial penalty curve. -// // Solved by sympy package: -// /* -// from sympy import * -// (x,a,b,c,d,r,z)=symbols('x a b c d r z') -// p = a + b*x + c*x*x + d*x*x*x -// p2 = p.subs(solve([p.subs(x, -r), p.diff(x).subs(x, -r), p.diff(x,x).subs(x, -r), p.subs(x, 0)-z], [a, b, c, d])) -// from sympy.plotting import plot -// plot(p2.subs(r,0.2).subs(z,1.), (x, -1, 3), adaptive=False, nb_of_points=400) -// */ -// if (overlap_distance < - nozzle_r) { -// // The extrudate is fully supported by the lower layer. This is the ideal case, therefore zero penalty. -// return 0.f; -// } else { -// float x = overlap_distance / nozzle_r; -// float x2 = x * x; -// float x3 = x2 * x; -// return weight_zero * (1.f + 3.f * x + 3.f * x2 + x3); -// } -// } -static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) +// ************ FOR BACKPORT COMPATIBILITY ONLY *************** +// Color mapping of a value into RGB false colors. +inline Vec3f value_to_rgbf(float minimum, float maximum, float value) { - return overlap_distance > nozzle_r ? weight_zero : 0.f; + 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 }; } - - -// 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) +// Color mapping of a value into RGB false colors. +inline Vec3i value_to_rgbi(float minimum, float maximum, float value) { - 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; + return (value_to_rgbf(minimum, maximum, value) * 255).cast(); +} +// *************************** + +template int sgn(T val) { + return int(T(0) < val) - int(val < T(0)); } +// base function: ((e^(((1)/(x^(2)+1)))-1)/(e-1)) +// checkout e.g. here: https://www.geogebra.org/calculator +float gauss(float value, float mean_x_coord, float mean_value, float falloff_speed) { + float shifted = value - mean_x_coord; + float denominator = falloff_speed * shifted * shifted + 1.0f; + float exponent = 1.0f / denominator; + return mean_value * (std::exp(exponent) - 1.0f) / (std::exp(1.0f) - 1.0f); +} +/// Coordinate frame +class Frame { +public: + Frame() { + mX = Vec3f(1, 0, 0); + mY = Vec3f(0, 1, 0); + mZ = Vec3f(0, 0, 1); + } -static Points::const_iterator project_point_to_polygon_and_insert(Polygon &polygon, const Point &pt, double eps) -{ - assert(polygon.points.size() >= 2); - if (polygon.points.size() <= 1) - if (polygon.points.size() == 1) - return polygon.points.begin(); + Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) : + mX(x), mY(y), mZ(z) { + } - Point pt_min; - double d_min = std::numeric_limits::max(); - size_t i_min = size_t(-1); + 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); + } - for (size_t i = 0; i < polygon.points.size(); ++ i) { - size_t j = i + 1; - if (j == polygon.points.size()) - j = 0; - const Point &p1 = polygon.points[i]; - const Point &p2 = polygon.points[j]; - const Slic3r::Point v_seg = p2 - p1; - const Slic3r::Point v_pt = pt - p1; - const int64_t l2_seg = int64_t(v_seg(0)) * int64_t(v_seg(0)) + int64_t(v_seg(1)) * int64_t(v_seg(1)); - int64_t t_pt = int64_t(v_seg(0)) * int64_t(v_pt(0)) + int64_t(v_seg(1)) * int64_t(v_pt(1)); - if (t_pt < 0) { - // Closest to p1. - double dabs = sqrt(int64_t(v_pt(0)) * int64_t(v_pt(0)) + int64_t(v_pt(1)) * int64_t(v_pt(1))); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - pt_min = p1; - } - } - else if (t_pt > l2_seg) { - // Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the next step. - continue; - } else { - // Closest to the segment. - assert(t_pt >= 0 && t_pt <= l2_seg); - int64_t d_seg = int64_t(v_seg(1)) * int64_t(v_pt(0)) - int64_t(v_seg(0)) * int64_t(v_pt(1)); - double d = double(d_seg) / sqrt(double(l2_seg)); - double dabs = std::abs(d); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - // Evaluate the foot point. - pt_min = p1; - double linv = double(d_seg) / double(l2_seg); - pt_min(0) = pt(0) - coord_t(floor(double(v_seg(1)) * linv + 0.5)); - pt_min(1) = pt(1) + coord_t(floor(double(v_seg(0)) * linv + 0.5)); - assert(Line(p1, p2).distance_to(pt_min) < scale_(1e-5)); - } + Vec3f to_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, size_t negative_volumes_start_index) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start"; + + //prepare uniform samples of a hemisphere + float step_size = 1.0f / SeamPlacer::sqr_rays_per_triangle; + std::vector precomputed_sample_directions( + SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); + for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_triangle; ++x_idx) { + float sample_x = x_idx * step_size + step_size / 2.0; + for (size_t y_idx = 0; y_idx < SeamPlacer::sqr_rays_per_triangle; ++y_idx) { + size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_triangle + y_idx; + float sample_y = y_idx * step_size + step_size / 2.0; + precomputed_sample_directions[dir_index] = sample_hemisphere_uniform( { sample_x, sample_y }); } } - assert(i_min != size_t(-1)); - if ((pt_min - polygon.points[i_min]).cast().norm() > eps) { - // Insert a new point on the segment i_min, i_min+1. - return polygon.points.insert(polygon.points.begin() + (i_min + 1), pt_min); - } - return polygon.points.begin() + i_min; -} + bool model_contains_negative_parts = negative_volumes_start_index < triangles.indices.size(); + std::vector result(triangles.indices.size()); + tbb::parallel_for(tbb::blocked_range(0, result.size()), + [&triangles, &precomputed_sample_directions, model_contains_negative_parts, negative_volumes_start_index, + &raycasting_tree, &result](tbb::blocked_range r) { + // Maintaining hits memory outside of the loop, so it does not have to be reallocated for each query. + std::vector hits; + for (size_t face_index = r.begin(); face_index < r.end(); ++face_index) { + FaceVisibilityInfo &dest = result[face_index]; + dest.visibility = 1.0f; + constexpr float decrease = 1.0f + / (SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); + Vec3i face = triangles.indices[face_index]; + Vec3f A = triangles.vertices[face.x()]; + Vec3f B = triangles.vertices[face.y()]; + Vec3f C = triangles.vertices[face.z()]; + Vec3f center = (A + B + C) / 3.0f; + Vec3f normal = ((B - A).cross(C - B)).normalized(); + // apply the local direction via Frame struct - the local_dir is with respect to +Z being forward + Frame f; + f.set_from_z(normal); -static std::vector polygon_angles_at_vertices(const Polygon &polygon, const std::vector &lengths, float min_arm_length) -{ - assert(polygon.points.size() + 1 == lengths.size()); - if (min_arm_length > 0.25f * lengths.back()) - min_arm_length = 0.25f * lengths.back(); - - // Find the initial prev / next point span. - size_t idx_prev = polygon.points.size(); - size_t idx_curr = 0; - size_t idx_next = 1; - while (idx_prev > idx_curr && lengths.back() - lengths[idx_prev] < min_arm_length) - -- idx_prev; - while (idx_next < idx_prev && lengths[idx_next] < min_arm_length) - ++ idx_next; - - std::vector angles(polygon.points.size(), 0.f); - for (; idx_curr < polygon.points.size(); ++ idx_curr) { - // Move idx_prev up until the distance between idx_prev and idx_curr is lower than min_arm_length. - if (idx_prev >= idx_curr) { - while (idx_prev < polygon.points.size() && lengths.back() - lengths[idx_prev] + lengths[idx_curr] > min_arm_length) - ++ idx_prev; - if (idx_prev == polygon.points.size()) - idx_prev = 0; - } - while (idx_prev < idx_curr && lengths[idx_curr] - lengths[idx_prev] > min_arm_length) - ++ idx_prev; - // Move idx_prev one step back. - if (idx_prev == 0) - idx_prev = polygon.points.size() - 1; - else - -- idx_prev; - // Move idx_next up until the distance between idx_curr and idx_next is greater than min_arm_length. - if (idx_curr <= idx_next) { - while (idx_next < polygon.points.size() && lengths[idx_next] - lengths[idx_curr] < min_arm_length) - ++ idx_next; - if (idx_next == polygon.points.size()) - idx_next = 0; - } - while (idx_next < idx_curr && lengths.back() - lengths[idx_curr] + lengths[idx_next] < min_arm_length) - ++ idx_next; - // Calculate angle between idx_prev, idx_curr, idx_next. - const Point &p0 = polygon.points[idx_prev]; - const Point &p1 = polygon.points[idx_curr]; - const Point &p2 = polygon.points[idx_next]; - const Point v1 = p1 - p0; - const Point v2 = p2 - p1; - int64_t dot = int64_t(v1(0))*int64_t(v2(0)) + int64_t(v1(1))*int64_t(v2(1)); - int64_t cross = int64_t(v1(0))*int64_t(v2(1)) - int64_t(v1(1))*int64_t(v2(0)); - float angle = float(atan2(double(cross), double(dot))); - angles[idx_curr] = angle; - } - - return angles; -} - - - -void SeamPlacer::init(const Print& print) -{ - m_enforcers.clear(); - m_blockers.clear(); - m_seam_history.clear(); - m_po_list.clear(); - - const std::vector& nozzle_dmrs = print.config().nozzle_diameter.values; - float max_nozzle_dmr = *std::max_element(nozzle_dmrs.begin(), nozzle_dmrs.end()); - - - std::vector temp_enf; - std::vector temp_blk; - std::vector temp_polygons; - - for (const PrintObject* po : print.objects()) { - - auto merge_and_offset = [po, &temp_polygons, max_nozzle_dmr](EnforcerBlockerType type, std::vector& out) { - // Offset the triangles out slightly. - auto offset_out = [](Polygon& input, float offset) -> ExPolygons { - ClipperLib::Paths out(1); - std::vector deltas(input.points.size(), offset); - input.make_counter_clockwise(); - out.front() = mittered_offset_path_scaled(input.points, deltas, 3.); - return ClipperPaths_to_Slic3rExPolygons(out, true); // perform union - }; - - - temp_polygons.clear(); - po->project_and_append_custom_facets(true, type, temp_polygons); - out.clear(); - out.reserve(temp_polygons.size()); - float offset = scale_(max_nozzle_dmr + po->config().elefant_foot_compensation); - for (Polygons &src : temp_polygons) { - out.emplace_back(ExPolygons()); - for (Polygon& plg : src) { - ExPolygons offset_explg = offset_out(plg, offset); - if (! offset_explg.empty()) - out.back().emplace_back(std::move(offset_explg.front())); - } - - offset = scale_(max_nozzle_dmr); - } - }; - merge_and_offset(EnforcerBlockerType::BLOCKER, temp_blk); - merge_and_offset(EnforcerBlockerType::ENFORCER, temp_enf); - - // Remember this PrintObject and initialize a store of enforcers and blockers for it. - m_po_list.push_back(po); - size_t po_idx = m_po_list.size() - 1; - m_enforcers.emplace_back(std::vector(temp_enf.size())); - m_blockers.emplace_back(std::vector(temp_blk.size())); - - // A helper class to store data to build the AABB tree from. - class CustomTriangleRef { - public: - CustomTriangleRef(size_t idx, - Point&& centroid, - BoundingBox&& bb) - : m_idx{idx}, m_centroid{centroid}, - m_bbox{AlignedBoxType(bb.min, bb.max)} - {} - size_t idx() const { return m_idx; } - const Point& centroid() const { return m_centroid; } - const TreeType::BoundingBox& bbox() const { return m_bbox; } - - private: - size_t m_idx; - Point m_centroid; - AlignedBoxType m_bbox; - }; - - // A lambda to extract the ExPolygons and save them into the member AABB tree. - // Will be called for enforcers and blockers separately. - auto add_custom = [](std::vector& src, std::vector& dest) { - // Go layer by layer, and append all the ExPolygons into the AABB tree. - size_t layer_idx = 0; - for (ExPolygons& expolys_on_layer : src) { - CustomTrianglesPerLayer& layer_data = dest[layer_idx]; - std::vector triangles_data; - layer_data.polys.reserve(expolys_on_layer.size()); - triangles_data.reserve(expolys_on_layer.size()); - - for (ExPolygon& expoly : expolys_on_layer) { - if (expoly.empty()) - continue; - layer_data.polys.emplace_back(std::move(expoly)); - triangles_data.emplace_back(layer_data.polys.size() - 1, - layer_data.polys.back().centroid(), - layer_data.polys.back().bounding_box()); - } - // All polygons are saved, build the AABB tree for them. - layer_data.tree.build(std::move(triangles_data)); - ++layer_idx; - } - }; - - add_custom(temp_enf, m_enforcers.at(po_idx)); - add_custom(temp_blk, m_blockers.at(po_idx)); - } -} - - - -void SeamPlacer::plan_perimeters(const std::vector perimeters, - const Layer& layer, SeamPosition seam_position, - Point last_pos, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid) -{ - // When printing the perimeters, we want the seams on external and internal perimeters to match. - // We have a list of perimeters in the order to be printed. Each internal perimeter must inherit - // the seam from the previous external perimeter. - - m_plan.clear(); - m_plan_idx = 0; - - if (perimeters.empty() || ! po) - return; - - m_plan.resize(perimeters.size()); - - for (int i = 0; i < int(perimeters.size()); ++i) { - if (perimeters[i]->role() == erExternalPerimeter && perimeters[i]->is_loop()) { - last_pos = this->calculate_seam( - layer, seam_position, *dynamic_cast(perimeters[i]), nozzle_dmr, - po, lower_layer_edge_grid, last_pos, false); - m_plan[i].external = true; - } - m_plan[i].seam_position = seam_position; - m_plan[i].layer = &layer; - m_plan[i].po = po; - m_plan[i].pt = last_pos; - } -} - - -void SeamPlacer::place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter, - const EdgeGrid::Grid* lower_layer_edge_grid) -{ - // const double seam_offset = nozzle_diameter; - - Point seam = last_pos; - if (! m_plan.empty() && m_plan_idx < m_plan.size()) { - if (m_plan[m_plan_idx].external) { - seam = m_plan[m_plan_idx].pt; - // One more heuristics: if the seam is too far from current nozzle position, - // try to place it again. This can happen in cases where the external perimeter - // does not belong to the preceding ones and they are ordered so they end up - // far from each other. - if ((seam.cast() - last_pos.cast()).squaredNorm() > std::pow(scale_(5.*nozzle_diameter), 2.)) - seam = this->calculate_seam(*m_plan[m_plan_idx].layer, m_plan[m_plan_idx].seam_position, loop, nozzle_diameter, - m_plan[m_plan_idx].po, lower_layer_edge_grid, last_pos, false); - - if (m_plan[m_plan_idx].seam_position == spAligned) - m_seam_history.add_seam(m_plan[m_plan_idx].po, m_plan[m_plan_idx].pt, loop.polygon().bounding_box()); - } - else { - if (!external_first) { - // Internal perimeter printed before the external. - // First get list of external seams. - std::vector ext_seams; - size_t external_cnt = 0; - for (size_t i = 0; i < m_plan.size(); ++i) { - if (m_plan[i].external) { - ext_seams.emplace_back(i); - ++external_cnt; - } - } - - if (!ext_seams.empty()) { - // First find the line segment closest to an external seam: - //int path_idx = 0; - //int line_idx = 0; - size_t ext_seam_idx = size_t(-1); - double min_dist_sqr = std::numeric_limits::max(); - std::vector lines_vect; - for (int i = 0; i < int(loop.paths.size()); ++i) { - lines_vect.emplace_back(loop.paths[i].polyline.lines()); - const Lines& lines = lines_vect.back(); - for (int j = 0; j < int(lines.size()); ++j) { - for (size_t k : ext_seams) { - double d_sqr = lines[j].distance_to_squared(m_plan[k].pt); - if (d_sqr < min_dist_sqr) { - //path_idx = i; - //line_idx = j; - ext_seam_idx = k; - min_dist_sqr = d_sqr; + for (const auto &dir : precomputed_sample_directions) { + Vec3f final_ray_dir = (f.to_world(dir)); + if (!model_contains_negative_parts) { + igl::Hit hitpoint; + // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and + // direction. + Vec3d final_ray_dir_d = final_ray_dir.cast(); + Vec3d ray_origin_d = (center + normal * 0.03).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) { + dest.visibility -= decrease; + } + } else { //TODO improve logic for order based boolean operations - consider order of volumes + Vec3d ray_origin_d = (center + normal * 0.1).cast(); // start above surface. + if (face_index >= negative_volumes_start_index) { // if casting from negative volume face, invert direction, change start pos + final_ray_dir = -1.0 * final_ray_dir; + ray_origin_d = (center - normal * 0.03).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) { + dest.visibility -= decrease; } } } } + } + }); - // Only accept seam that is reasonably close. - if (ext_seam_idx != size_t(-1)) { - // How many nozzle diameters is considered "close"? - const double nozzle_d_limit = 2. * (1. + m_plan.size() / external_cnt); - const double limit_dist_sqr = double(scale_(scale_((unscale(m_plan[ext_seam_idx].pt) - unscale(m_plan[m_plan_idx].pt)).squaredNorm() * std::pow(nozzle_d_limit * nozzle_diameter, 2.)))); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end"; - if (min_dist_sqr < limit_dist_sqr) { - // Now find a projection of the external seam - //const Lines& lines = lines_vect[path_idx]; - //Point closest = m_plan[ext_seam_idx].pt.projection_onto(lines[line_idx]); + return result; +} - // This code does staggering of internal perimeters, turned off for now. - // - // double dist = (closest.cast() - lines[line_idx].b.cast()).norm(); - // - // // And walk along the perimeter until we make enough space for - // // seams of all perimeters beforethe external one. - // double offset = (ext_seam_idx - m_plan_idx) * scale_(seam_offset); - // double last_offset = offset; - // offset -= dist; - // const Point* a = &closest; - // const Point* b = &lines[line_idx].b; - // while (++line_idx < int(lines.size()) && offset > 0.) { - // last_offset = offset; - // offset -= lines[line_idx].length(); - // a = &lines[line_idx].a; - // b = &lines[line_idx].b; - // } - // - // // We have walked far enough, too far maybe. Interpolate on the - // // last segment to find the end precisely. - // offset = std::min(0., offset); // In case that offset is still positive (we may have "wrapped around") - // double ratio = last_offset / (last_offset - offset); - // seam = (a->cast() + ((b->cast() - a->cast()) * ratio)).cast(); - seam = m_plan[ext_seam_idx].pt; - } +std::vector calculate_polygon_angles_at_vertices(const Polygon &polygon, const std::vector &lengths, + float min_arm_length) { + std::vector result(polygon.size()); + + if (polygon.size() == 1) { + result[0] = 0.0f; + } + + size_t idx_prev = 0; + size_t idx_curr = 0; + size_t idx_next = 0; + + float distance_to_prev = 0; + float distance_to_next = 0; + + //push idx_prev far enough back as initialization + while (distance_to_prev < min_arm_length) { + idx_prev = Slic3r::prev_idx_modulo(idx_prev, polygon.size()); + distance_to_prev += lengths[idx_prev]; + } + + for (size_t _i = 0; _i < polygon.size(); ++_i) { + // pull idx_prev to current as much as possible, while respecting the min_arm_length + while (distance_to_prev - lengths[idx_prev] > min_arm_length) { + distance_to_prev -= lengths[idx_prev]; + idx_prev = Slic3r::next_idx_modulo(idx_prev, polygon.size()); + } + + //push idx_next forward as far as needed + while (distance_to_next < min_arm_length) { + distance_to_next += lengths[idx_next]; + idx_next = Slic3r::next_idx_modulo(idx_next, polygon.size()); + } + + // Calculate angle between idx_prev, idx_curr, idx_next. + const Point &p0 = polygon.points[idx_prev]; + const Point &p1 = polygon.points[idx_curr]; + const Point &p2 = polygon.points[idx_next]; + result[idx_curr] = float(angle(p1 - p0, p2 - p1)); + + // increase idx_curr by one + float curr_distance = lengths[idx_curr]; + idx_curr++; + distance_to_prev += curr_distance; + distance_to_next -= curr_distance; + } + + return result; +} + +// structure to store global information about the model - occlusion hits, enforcers, blockers +struct GlobalModelInfo { + indexed_triangle_set model; + std::vector triangle_neighbours; + AABBTreeIndirect::Tree<3, float> model_tree; + std::vector visiblity_info; + indexed_triangle_set enforcers; + indexed_triangle_set blockers; + AABBTreeIndirect::Tree<3, float> enforcers_tree; + AABBTreeIndirect::Tree<3, float> blockers_tree; + + bool is_enforced(const Vec3f &position, float radius) const { + if (enforcers.empty()) { + return false; + } + float radius_sqr = radius * radius; + return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices, + enforcers_tree, position, radius_sqr); + } + + bool is_blocked(const Vec3f &position, float radius) const { + if (blockers.empty()) { + return false; + } + float radius_sqr = radius * radius; + return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices, + blockers_tree, position, radius_sqr); + } + + float calculate_point_visibility(const Vec3f &position) const { + size_t hit_idx; + Vec3f hit_point; + if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree, + position, hit_idx, hit_point) >= 0) { + float visibility = visiblity_info[hit_idx].visibility; + Vec3i neighbours = this->triangle_neighbours[hit_idx]; + size_t n_count = 0; + for (int neighbour : neighbours) { + if (neighbour >= 0) { + visibility += visiblity_info[neighbour].visibility; + n_count++; + } + } + return visibility / (1 + n_count); + } else { + return 0.0f; + } + } + +#ifdef DEBUG_FILES + void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const { + indexed_triangle_set divided_mesh = obj_mesh; + Slic3r::CNumericLocalesSetter locales_setter; + + FILE *fp = boost::nowide::fopen(file_name, "w"); + if (fp == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << file_name << " for writing"; + return; + } + + for (size_t i = 0; i < divided_mesh.vertices.size(); ++i) { + float visibility = calculate_point_visibility(divided_mesh.vertices[i]); + Vec3f color = value_to_rgbf(0.0f, 1.0f, visibility); + fprintf(fp, "v %f %f %f %f %f %f\n", + divided_mesh.vertices[i](0), divided_mesh.vertices[i](1), divided_mesh.vertices[i](2), + color(0), color(1), color(2)); + } + for (size_t i = 0; i < divided_mesh.indices.size(); ++i) + fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1, + divided_mesh.indices[i][2] + 1); + fclose(fp); + + } +#endif +} +; + +//Extract perimeter polygons of the given layer +Polygons extract_perimeter_polygons(const Layer *layer, const SeamPosition configured_seam_preference, + std::vector &corresponding_regions_out) { + Polygons polygons; + for (const LayerRegion *layer_region : layer->regions()) { + for (const ExtrusionEntity *ex_entity : layer_region->perimeters.entities) { + if (ex_entity->is_collection()) { //collection of inner, outer, and overhang perimeters + for (const ExtrusionEntity *perimeter : static_cast(ex_entity)->entities) { + if (perimeter->role() == ExtrusionRole::erExternalPerimeter + || (perimeter->role() == ExtrusionRole::erPerimeter + && configured_seam_preference == spRandom)) { //for random seam alignment, extract all perimeters + Points p; + perimeter->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); } } + if (polygons.empty()) { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); + } + } else { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(std::move(p)); + corresponding_regions_out.push_back(layer_region); } - else { - // We should have a candidate ready from before. If not, use last_pos. - if (m_plan_idx > 0 && m_plan[m_plan_idx - 1].precalculated) - seam = m_plan[m_plan_idx - 1].pt; - } - - // seam now contains a hot candidate for internal seam. Use it unless there is a sharp corner nearby. - // We will call the normal seam planning function, pretending that we are currently at the candidate point - // and set to spNearest. If the ideal seam it finds is close to current candidate, use it. - // This is to prevent having seams very close to corners, just because of external seam position. - seam = calculate_seam(*m_plan[m_plan_idx].layer, spNearest, loop, nozzle_diameter, - m_plan[m_plan_idx].po, lower_layer_edge_grid, seam, true); } - m_plan[m_plan_idx].pt = seam; } - - // Split the loop at the point with a minium penalty. - if (!loop.split_at_vertex(seam)) - // The point is not in the original loop. Insert it. - loop.split_at(seam, true); - - if (external_first && m_plan_idx+1 1) { -// const ExtrusionPath& last = loop.paths.back(); -// auto it = last.polyline.points.crbegin() + 1; -// for (; it != last.polyline.points.crend(); ++it) { -// running_sqr += (it->cast() - (it - 1)->cast()).squaredNorm(); -// if (running_sqr > dist_sqr) -// break; -// running_sqr_last = running_sqr; -// } -// if (running_sqr <= dist_sqr) -// it = last.polyline.points.crend() - 1; -// // Now interpolate. -// double ratio = (std::sqrt(dist_sqr) - std::sqrt(running_sqr_last)) / (std::sqrt(running_sqr) - std::sqrt(running_sqr_last)); -// m_plan[m_plan_idx + 1].pt = ((it - 1)->cast() + (it->cast() - (it - 1)->cast()) * std::min(ratio, 1.)).cast(); -// m_plan[m_plan_idx + 1].precalculated = true; - m_plan[m_plan_idx + 1].pt = m_plan[m_plan_idx].pt; - m_plan[m_plan_idx + 1].precalculated = true; -// } + if (polygons.empty()) { // If there are no perimeter polygons for whatever reason (disabled perimeters .. ) insert dummy point + // it is easier than checking everywhere if the layer is not emtpy, no seam will be placed to this layer anyway + polygons.emplace_back(std::vector { Point { 0, 0 } }); + corresponding_regions_out.push_back(nullptr); } - ++m_plan_idx; + return polygons; } +// Insert SeamCandidates created from perimeter polygons in to the result vector. +// Compute its type (Enfrocer,Blocker), angle, and position +//each SeamCandidate also contains pointer to shared Perimeter structure representing the polygon +// if Custom Seam modifiers are present, oversamples the polygon if necessary to better fit user intentions +void process_perimeter_polygon(const Polygon &orig_polygon, float z_coord, const LayerRegion *region, + const GlobalModelInfo &global_model_info, PrintObjectSeamData::LayerSeams &result) { + if (orig_polygon.size() == 0) { + return; + } -// Returns "best" seam for a given perimeter. -Point SeamPlacer::calculate_seam(const Layer& layer, const SeamPosition seam_position, - const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest) -{ - Polygon polygon = loop.polygon(); + Polygon polygon = orig_polygon; bool was_clockwise = polygon.make_counter_clockwise(); - const coord_t nozzle_r = coord_t(scale_(0.5 * nozzle_dmr) + 0.5); - size_t po_idx = std::find(m_po_list.begin(), m_po_list.end(), po) - m_po_list.begin(); - - // Find current layer in respective PrintObject. Cache the result so the - // lookup is only done once per layer, not for each loop. - const Layer* layer_po = nullptr; - if (po == m_last_po && layer.print_z == m_last_print_z) - layer_po = m_last_layer_po; - else { - layer_po = po ? po->get_layer_at_printz(layer.print_z) : nullptr; - m_last_po = po; - m_last_print_z = layer.print_z; - m_last_layer_po = layer_po; + std::vector lengths { }; + for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) { + lengths.push_back(std::max((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm(), 0.01)); } - if (! layer_po) - return last_pos; + lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.01)); - // Index of this layer in the respective PrintObject. - size_t layer_idx = layer_po->id() - po->layers().front()->id(); // raft layers + std::vector local_angles = calculate_polygon_angles_at_vertices(polygon, lengths, + SeamPlacer::polygon_local_angles_arm_distance); - assert(layer_idx < po->layer_count()); + result.perimeters.push_back( { }); + Perimeter &perimeter = result.perimeters.back(); - const bool custom_seam = loop.role() == erExternalPerimeter && this->is_custom_seam_on_layer(layer_idx, po_idx); - - if (custom_seam) { - // Seam enf/blockers can begin and end in between the original vertices. - // Let add extra points in between and update the leghths. - polygon.densify(MINIMAL_POLYGON_SIDE); + std::queue orig_polygon_points { }; + for (size_t index = 0; index < polygon.size(); ++index) { + Vec2f unscaled_p = unscale(polygon[index]).cast(); + orig_polygon_points.emplace(unscaled_p.x(), unscaled_p.y(), z_coord); } + Vec3f first = orig_polygon_points.front(); + std::queue oversampled_points { }; + size_t orig_angle_index = 0; + perimeter.start_index = result.points.size(); + perimeter.flow_width = region != nullptr ? region->flow(FlowRole::frExternalPerimeter).width() : 0.0f; + bool some_point_enforced = false; + while (!orig_polygon_points.empty() || !oversampled_points.empty()) { + EnforcedBlockedSeamPoint type = EnforcedBlockedSeamPoint::Neutral; + Vec3f position; + float local_ccw_angle = 0; + bool orig_point = false; + if (!oversampled_points.empty()) { + position = oversampled_points.front(); + oversampled_points.pop(); + } else { + position = orig_polygon_points.front(); + orig_polygon_points.pop(); + local_ccw_angle = was_clockwise ? -local_angles[orig_angle_index] : local_angles[orig_angle_index]; + orig_angle_index++; + orig_point = true; + } - if (seam_position != spRandom) { - // Retrieve the last start position for this object. - float last_pos_weight = 1.f; + if (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { + type = EnforcedBlockedSeamPoint::Enforced; + some_point_enforced = true; + } - if (seam_position == spAligned) { - // Seam is aligned to the seam at the preceding layer. - if (po != nullptr) { - std::optional pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, loop.polygon().bounding_box()); - if (pos.has_value()) { - last_pos = *pos; - last_pos_weight = is_custom_enforcer_on_layer(layer_idx, po_idx) ? 0.f : 1.f; + if (global_model_info.is_blocked(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { + type = EnforcedBlockedSeamPoint::Blocked; + } + + if (orig_point) { + Vec3f pos_of_next = orig_polygon_points.empty() ? first : orig_polygon_points.front(); + float distance_to_next = (position - pos_of_next).norm(); + if (global_model_info.is_enforced(position, distance_to_next)) { + Vec3f vec_to_next = (pos_of_next - position).normalized(); + float step_size = SeamPlacer::enforcer_oversampling_distance; + float step = step_size; + while (step < distance_to_next) { + oversampled_points.push(position + vec_to_next * step); + step += step_size; } } } - else if (seam_position == spRear) { - // Object is centered around (0,0) in its current coordinate system. - last_pos.x() = 0; - last_pos.y() = coord_t(3. * po->bounding_box().radius()); - last_pos_weight = 5.f; - } if (seam_position == spNearest) { - // last_pos already contains current nozzle position - } - // Insert a projection of last_pos into the polygon. - size_t last_pos_proj_idx; - { - auto it = project_point_to_polygon_and_insert(polygon, last_pos, 0.1 * nozzle_r); - last_pos_proj_idx = it - polygon.points.begin(); - } - - // Parametrize the polygon by its length. - std::vector lengths = polygon.parameter_by_length(); - - // For each polygon point, store a penalty. - // First calculate the angles, store them as penalties. The angles are caluculated over a minimum arm length of nozzle_r. - std::vector penalties = polygon_angles_at_vertices(polygon, lengths, float(nozzle_r)); - // No penalty for reflex points, slight penalty for convex points, high penalty for flat surfaces. - const float penaltyConvexVertex = 1.f; - const float penaltyFlatSurface = 3.f; - const float penaltyOverhangHalf = 10.f; - // Penalty for visible seams. - for (size_t i = 0; i < polygon.points.size(); ++ i) { - float ccwAngle = penalties[i]; - if (was_clockwise) - ccwAngle = - ccwAngle; - float penalty = 0; - if (ccwAngle <- float(0.6 * PI)) - // Sharp reflex vertex. We love that, it hides the seam perfectly. - penalty = 0.f; - else if (ccwAngle > float(0.6 * PI)) - // Seams on sharp convex vertices are more visible than on reflex vertices. - penalty = penaltyConvexVertex; - else if (ccwAngle < 0.f) { - // Interpolate penalty between maximum and zero. - penalty = penaltyFlatSurface * bspline_kernel(ccwAngle * float(PI * 2. / 3.)); - } else { - assert(ccwAngle >= 0.f); - // Interpolate penalty between maximum and the penalty for a convex vertex. - penalty = penaltyConvexVertex + (penaltyFlatSurface - penaltyConvexVertex) * bspline_kernel(ccwAngle * float(PI * 2. / 3.)); - } - // Give a negative penalty for points close to the last point or the prefered seam location. - float dist_to_last_pos_proj = (i < last_pos_proj_idx) ? - std::min(lengths[last_pos_proj_idx] - lengths[i], lengths.back() - lengths[last_pos_proj_idx] + lengths[i]) : - std::min(lengths[i] - lengths[last_pos_proj_idx], lengths.back() - lengths[i] + lengths[last_pos_proj_idx]); - float dist_max = 0.1f * lengths.back(); // 5.f * nozzle_dmr - penalty -= last_pos_weight * bspline_kernel(dist_to_last_pos_proj / dist_max); - penalties[i] = std::max(0.f, penalty); - if (prefer_nearest) { - // This hack limits the search around the nearest position projection. - penalties[i] += dist_to_last_pos_proj > 6.f * nozzle_r ? 100.f : 0.f; - } - } - - // Penalty for overhangs. - if (lower_layer_edge_grid) { - // Use the edge grid distance field structure over the lower layer to calculate overhangs. - coord_t nozzle_r = coord_t(std::floor(scale_(0.5 * nozzle_dmr) + 0.5)); - coord_t search_r = coord_t(std::floor(scale_(0.8 * nozzle_dmr) + 0.5)); - for (size_t i = 0; i < polygon.points.size(); ++ i) { - const Point &p = polygon.points[i]; - coordf_t dist; - // Signed distance is positive outside the object, negative inside the object. - // The point is considered at an overhang, if it is more than nozzle radius - // outside of the lower layer contour. - [[maybe_unused]] bool found = lower_layer_edge_grid->signed_distance(p, search_r, dist); - // If the approximate Signed Distance Field was initialized over lower_layer_edge_grid, - // then the signed distnace shall always be known. - assert(found); - penalties[i] += extrudate_overlap_penalty(float(nozzle_r), penaltyOverhangHalf, float(dist)); - } - } - - // Custom seam. Huge (negative) constant penalty is applied inside - // blockers (enforcers) to rule out points that should not win. - if (custom_seam) - this->apply_custom_seam(polygon, po_idx, penalties, lengths, layer_idx, seam_position); - - // Find a point with a minimum penalty. - size_t idx_min = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); - - if (seam_position != spAligned || ! is_custom_enforcer_on_layer(layer_idx, po_idx)) { - // Very likely the weight of idx_min is very close to the weight of last_pos_proj_idx. - // In that case use last_pos_proj_idx instead. - float penalty_aligned = penalties[last_pos_proj_idx]; - float penalty_min = penalties[idx_min]; - float penalty_diff_abs = std::abs(penalty_min - penalty_aligned); - float penalty_max = std::max(std::abs(penalty_min), std::abs(penalty_aligned)); - float penalty_diff_rel = (penalty_max == 0.f) ? 0.f : penalty_diff_abs / penalty_max; - // printf("Align seams, penalty aligned: %f, min: %f, diff abs: %f, diff rel: %f\n", penalty_aligned, penalty_min, penalty_diff_abs, penalty_diff_rel); - if (std::abs(penalty_diff_rel) < 0.05) { - // Penalty of the aligned point is very close to the minimum penalty. - // Align the seams as accurately as possible. - idx_min = last_pos_proj_idx; - } - } - - - // Export the contour into a SVG file. - #if 0 - { - static int iRun = 0; - SVG svg(debug_out_path("GCode_extrude_loop-%d.svg", iRun ++)); - if (m_layer->lower_layer != NULL) - svg.draw(m_layer->lower_layer->slices); - for (size_t i = 0; i < loop.paths.size(); ++ i) - svg.draw(loop.paths[i].as_polyline(), "red"); - Polylines polylines; - for (size_t i = 0; i < loop.paths.size(); ++ i) - polylines.push_back(loop.paths[i].as_polyline()); - Slic3r::Polygons polygons; - coordf_t nozzle_dmr = EXTRUDER_CONFIG(nozzle_diameter); - coord_t delta = scale_(0.5*nozzle_dmr); - Slic3r::offset(polylines, &polygons, delta); -// for (size_t i = 0; i < polygons.size(); ++ i) svg.draw((Polyline)polygons[i], "blue"); - svg.draw(last_pos, "green", 3); - svg.draw(polygon.points[idx_min], "yellow", 3); - svg.Close(); - } - #endif - return polygon.points[idx_min]; - - } else - return this->get_random_seam(layer_idx, polygon, po_idx); -} - - -Point SeamPlacer::get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_idx, - bool* saw_custom) const -{ - // Parametrize the polygon by its length. - const std::vector lengths = polygon.parameter_by_length(); - - // Which of the points are inside enforcers/blockers? - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_idx, polygon, po_idx, enforcers_idxs, blockers_idxs); - - bool has_enforcers = ! enforcers_idxs.empty(); - bool has_blockers = ! blockers_idxs.empty(); - if (saw_custom) - *saw_custom = has_enforcers || has_blockers; - - assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end())); - assert(std::is_sorted(blockers_idxs.begin(), blockers_idxs.end())); - std::vector edges; - - // Lambda to calculate lengths of all edges of interest. Last parameter - // decides whether to measure edges inside or outside idxs. - // Negative number = not an edge of interest. - auto get_valid_length = [&lengths](const std::vector& idxs, - std::vector& edges, - bool measure_inside_edges) -> float - { - // First mark edges we are interested in by assigning a positive number. - edges.assign(lengths.size()-1, measure_inside_edges ? -1.f : 1.f); - for (size_t i=0; i the edge between them is the enforcer/blocker. - bool inside_edge = ((i != idxs.size()-1 && idxs[i+1] == this_pt_idx + 1) - || (i == idxs.size()-1 && idxs.back() == lengths.size()-2 && idxs[0] == 0)); - if (inside_edge) - edges[this_pt_idx] = measure_inside_edges ? 1.f : -1.f; - } - // Now measure them. - float running_total = 0.f; - for (size_t i=0; i 0.f) { - edges[i] = lengths[i+1] - lengths[i]; - running_total += edges[i]; - } - } - return running_total; - }; - - // Find all seam candidate edges and their lengths. - float valid_length = 0.f; - if (has_enforcers) - valid_length = get_valid_length(enforcers_idxs, edges, true); - - if (! has_enforcers || valid_length == 0.f) { - // Second condition covers case with isolated enf points. Given how the painted - // triangles are projected, this should not happen. Stay on the safe side though. - if (has_blockers) - valid_length = get_valid_length(blockers_idxs, edges, false); - if (valid_length == 0.f) // No blockers or everything blocked - use the whole polygon. - valid_length = lengths.back(); + result.points.emplace_back(position, perimeter, local_ccw_angle, type); } - assert(valid_length != 0.f); - // Now generate a random length and find the respective edge. - float rand_len = valid_length * (rand()/float(RAND_MAX)); - size_t pt_idx = 0; // Index of the edge where to put the seam. - if (valid_length == lengths.back()) { - // Whole polygon is used for placing the seam. - auto it = std::lower_bound(lengths.begin(), lengths.end(), rand_len); - pt_idx = it == lengths.begin() ? 0 : (it-lengths.begin()-1); // this takes care of a corner case where rand() returns 0 - } else { - float running = 0.f; - for (size_t i=0; i 0.f ? edges[i] : 0.f; - if (running >= rand_len) { - pt_idx = i; + + perimeter.end_index = result.points.size() - 1; + + // We will find first patch of enforced points (patch: continuous section of enforced points) and select the middle + // point, which will have priority during alignment + // If there are multiple enforced patches in the perimeter, others are ignored + if (some_point_enforced) { + size_t perimeter_size = perimeter.end_index - perimeter.start_index + 1; + const auto next_index = [&](size_t idx) { + return perimeter.start_index + Slic3r::next_idx_modulo(idx - perimeter.start_index, perimeter_size); + }; + + size_t first_enforced_idx = perimeter.start_index; + for (size_t _ = 0; _ < perimeter_size; ++_) { + if (result.points[first_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced && + result.points[next_index(first_enforced_idx)].type == EnforcedBlockedSeamPoint::Enforced) { break; } + first_enforced_idx = next_index(first_enforced_idx); } - } + first_enforced_idx = next_index(first_enforced_idx); - if (! has_enforcers && ! has_blockers) { - // The polygon may be too coarse, calculate the point exactly. - assert(valid_length == lengths.back()); - bool last_seg = pt_idx == polygon.points.size()-1; - size_t next_idx = last_seg ? 0 : pt_idx+1; - const Point& prev = polygon.points[pt_idx]; - const Point& next = polygon.points[next_idx]; - assert(next_idx == 0 || pt_idx+1 == next_idx); - coordf_t diff_x = next.x() - prev.x(); - coordf_t diff_y = next.y() - prev.y(); - coordf_t dist = lengths[last_seg ? pt_idx+1 : next_idx] - lengths[pt_idx]; - return Point(prev.x() + (rand_len - lengths[pt_idx]) * (diff_x/dist), - prev.y() + (rand_len - lengths[pt_idx]) * (diff_y/dist)); - - } else { - // The polygon should be dense enough. - return polygon.points[pt_idx]; - } -} - - - - - - - - -void SeamPlacer::get_enforcers_and_blockers(size_t layer_id, - const Polygon& polygon, - size_t po_idx, - std::vector& enforcers_idxs, - std::vector& blockers_idxs) const -{ - enforcers_idxs.clear(); - blockers_idxs.clear(); - - auto is_inside = [](const Point& pt, - const CustomTrianglesPerLayer& custom_data) -> bool { - assert(! custom_data.polys.empty()); - // Now ask the AABB tree which polygons we should check and check them. - std::vector candidates; - AABBTreeIndirect::get_candidate_idxs(custom_data.tree, pt, candidates); - if (! candidates.empty()) - for (size_t idx : candidates) - if (custom_data.polys[idx].contains(pt)) - return true; - return false; - }; - - if (! m_enforcers[po_idx].empty()) { - const CustomTrianglesPerLayer& enforcers = m_enforcers[po_idx][layer_id]; - if (! enforcers.polys.empty()) { - for (size_t i=0; i orig_large_angle_points_indices { }; + std::vector viable_points_indices { }; + size_t last_enforced_idx = first_enforced_idx; + for (size_t _ = 0; _ < perimeter_size; ++_) { + if (result.points[last_enforced_idx].type != EnforcedBlockedSeamPoint::Enforced) { + break; } + viable_points_indices.push_back(last_enforced_idx); + if (abs(result.points[last_enforced_idx].local_ccw_angle) > 0.3 * PI) { + orig_large_angle_points_indices.push_back(last_enforced_idx); + } + last_enforced_idx = next_index(last_enforced_idx); } - } - - if (! m_blockers[po_idx].empty()) { - const CustomTrianglesPerLayer& blockers = m_blockers[po_idx][layer_id]; - if (! blockers.polys.empty()) { - for (size_t i=0; i find_enforcer_centers(const Polygon& polygon, - const std::vector& lengths, - const std::vector& enforcers_idxs) -{ - std::vector out; - assert(polygon.points.size()+1 == lengths.size()); - assert(std::is_sorted(enforcers_idxs.begin(), enforcers_idxs.end())); - if (polygon.size() < 2 || enforcers_idxs.empty()) - return out; - - auto get_center_idx = [&lengths](size_t start_idx, size_t end_idx) -> size_t { - assert(end_idx >= start_idx); - if (start_idx == end_idx) - return start_idx; - float t_c = lengths[start_idx] + 0.5f * (lengths[end_idx] - lengths[start_idx]); - auto it = std::lower_bound(lengths.begin() + start_idx, lengths.begin() + end_idx, t_c); - int ret = it - lengths.begin(); - return ret; - }; - - int last_enforcer_start_idx = enforcers_idxs.front(); - bool first_pt_in_list = enforcers_idxs.front() != 0; - bool last_pt_in_list = enforcers_idxs.back() == polygon.points.size() - 1; - bool wrap_around = last_pt_in_list && first_pt_in_list; - - for (size_t i=0; i 0); + if (orig_large_angle_points_indices.empty()) { + size_t central_idx = viable_points_indices[viable_points_indices.size() / 2]; + result.points[central_idx].central_enforcer = true; } else { - if (! wrap_around) { - // we can safely use the last enforcer point. - out.push_back(get_center_idx(last_enforcer_start_idx, enforcers_idxs[i])); + size_t central_idx = orig_large_angle_points_indices.size() / 2; + result.points[orig_large_angle_points_indices[central_idx]].central_enforcer = true; + } + } + +} + +// Get index of previous and next perimeter point of the layer. Because SeamCandidates of all polygons of the given layer +// are sequentially stored in the vector, each perimeter contains info about start and end index. These vales are used to +// deduce index of previous and next neigbour in the corresponding perimeter. +std::pair find_previous_and_next_perimeter_point(const std::vector &perimeter_points, + size_t point_index) { + const SeamCandidate ¤t = perimeter_points[point_index]; + int prev = point_index - 1; //for majority of points, it is true that neighbours lie behind and in front of them in the vector + int next = point_index + 1; + + if (point_index == current.perimeter.start_index) { + // if point_index is equal to start, it means that the previous neighbour is at the end + prev = current.perimeter.end_index; + } + + if (point_index == current.perimeter.end_index) { + // if point_index is equal to end, than next neighbour is at the start + next = current.perimeter.start_index; + } + + assert(prev >= 0); + assert(next >= 0); + return {size_t(prev),size_t(next)}; +} + +// Computes all global model info - transforms object, performs raycasting +void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather occlusion meshes: start"; + auto obj_transform = po->trafo_centered(); + indexed_triangle_set triangle_set; + indexed_triangle_set negative_volumes_set; + //add all parts + for (const ModelVolume *model_volume : po->model_object()->volumes) { + if (model_volume->type() == ModelVolumeType::MODEL_PART + || model_volume->type() == ModelVolumeType::NEGATIVE_VOLUME) { + auto model_transformation = model_volume->get_matrix(); + indexed_triangle_set model_its = model_volume->mesh().its; + its_transform(model_its, model_transformation); + if (model_volume->type() == ModelVolumeType::MODEL_PART) { + its_merge(triangle_set, model_its); + } else { + its_merge(negative_volumes_set, model_its); } } } + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather occlusion meshes: end"; - if (wrap_around) { - // Update first center already found. - if (out.empty()) { - // Probably an enforcer around the whole contour. Return nothing. - return out; - } + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: simplify occlusion meshes: start"; - // find last point of the enforcer at the beginning: - size_t idx = 0; - while (enforcers_idxs[idx]+1 == enforcers_idxs[idx+1]) - ++idx; + //simplify raycasting mesh + { + its_quadric_edge_collapse(triangle_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr, + nullptr, + nullptr); + float triangle_set_area = tbb::parallel_reduce(tbb::blocked_range(0, triangle_set.indices.size()), 0, + [&triangle_set]( + tbb::blocked_range r, float sum) { + 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()]; + sum += 0.5f * (b - a).cross(c - a).norm(); + } + return sum; + }, std::plus()); - float t_s = lengths[last_enforcer_start_idx]; - float t_e = lengths[idx]; - float half_dist = 0.5f * (t_e + lengths.back() - t_s); - float t_c = (half_dist > t_e) ? t_s + half_dist : t_e - half_dist; - - 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); + float target_triangle_area = triangle_set_area / SeamPlacer::raycasting_subdivision_target_triangle_count; + float target_triangle_length = 2 * 1.316 * sqrtf(target_triangle_area); //assuming 30-30-120 triangle + float subdivision_length = std::max(SeamPlacer::raycasting_subdivision_target_length, target_triangle_length); + triangle_set = its_subdivide(triangle_set, subdivision_length); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: triangle set after subdivision: " << triangle_set.indices.size(); } - return out; + + //simplify negative volumes + { + its_quadric_edge_collapse(negative_volumes_set, SeamPlacer::raycasting_decimation_target_triangle_count, + nullptr, + nullptr, + nullptr); + float negative_volumes_set_area = tbb::parallel_reduce( + tbb::blocked_range(0, negative_volumes_set.indices.size()), 0, + [&negative_volumes_set]( + tbb::blocked_range r, float sum) { + for (size_t t_idx = r.begin(); t_idx < r.end(); ++t_idx) { + const Vec3f &a = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].x()]; + const Vec3f &b = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].y()]; + const Vec3f &c = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].z()]; + sum += 0.5f * (b - a).cross(c - a).norm(); + } + return sum; + }, std::plus()); + + float target_triangle_area = negative_volumes_set_area + / SeamPlacer::raycasting_subdivision_target_triangle_count; + float target_triangle_length = 2 * 1.316 * sqrtf(target_triangle_area); //assuming 30-30-120 triangle + float subdivision_length = std::max(SeamPlacer::raycasting_subdivision_target_length, target_triangle_length); + negative_volumes_set = its_subdivide(negative_volumes_set, subdivision_length); + } + + size_t negative_volumes_start_index = triangle_set.indices.size(); + its_merge(triangle_set, negative_volumes_set); + its_transform(triangle_set, obj_transform); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: simplify occlusion meshes: end"; + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer:build AABB tree: start"; + auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices, + triangle_set.indices); + + std::vector neighbours = its_face_neighbors_par(triangle_set); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer:build AABB tree: end"; + result.model = triangle_set; + result.triangle_neighbours = neighbours; + result.model_tree = raycasting_tree; + result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set, negative_volumes_start_index); + +#ifdef DEBUG_FILES + auto filename = debug_out_path("visiblity.obj"); + result.debug_export(triangle_set, filename.c_str()); +#endif } +void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: start"; + auto obj_transform = po->trafo_centered(); -void SeamPlacer::apply_custom_seam(const Polygon& polygon, size_t po_idx, - std::vector& penalties, - const std::vector& lengths, - int layer_id, SeamPosition seam_position) const -{ - if (! is_custom_seam_on_layer(layer_id, po_idx)) - return; + for (const ModelVolume *mv : po->model_object()->volumes) { + if (mv->is_seam_painted()) { + auto model_transformation = obj_transform * mv->get_matrix(); - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs); + indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER); + its_transform(enforcers, model_transformation); + its_merge(result.enforcers, enforcers); - for (size_t i : enforcers_idxs) { - assert(i < penalties.size()); - penalties[i] -= float(ENFORCER_BLOCKER_PENALTY); - } - for (size_t i : blockers_idxs) { - assert(i < penalties.size()); - penalties[i] += float(ENFORCER_BLOCKER_PENALTY); - } - if (seam_position == spAligned) { - std::vector enf_centers = find_enforcer_centers(polygon, lengths, enforcers_idxs); - for (size_t idx : enf_centers) { - assert(idx < penalties.size()); - penalties[idx] += ENFORCER_CENTER_PENALTY; + indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER); + its_transform(blockers, model_transformation); + its_merge(result.blockers, blockers); } } -#if 0 - std::ostringstream os; - os << std::setw(3) << std::setfill('0') << layer_id; - int a = scale_(30.); - SVG svg("custom_seam" + os.str() + ".svg", BoundingBox(Point(-a, -a), Point(a, a))); - if (! m_enforcers[po_idx].empty()) - svg.draw(m_enforcers[po_idx][layer_id].polys, "blue"); - if (! m_blockers[po_idx].empty()) - svg.draw(m_blockers[po_idx][layer_id].polys, "red"); + result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices, + result.enforcers.indices); + result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices, + result.blockers.indices); - if (! blockers_idxs.empty()) { - ; + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end"; +} + +struct SeamComparator { + SeamPosition setup; + SeamComparator(SeamPosition setup) : + setup(setup) { } + // Standard comparator, must respect the requirements of comparators (e.g. give same result on same inputs) for sorting usage + // should return if a is better seamCandidate than b + bool is_first_better(const SeamCandidate &a, const SeamCandidate &b, const Vec2f &preffered_location = Vec2f { 0.0f, + 0.0f }) const { + if (setup == SeamPosition::spAligned && a.central_enforcer != b.central_enforcer) { + return a.central_enforcer; + } - size_t min_idx = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); + // Blockers/Enforcers discrimination, top priority + if (a.type != b.type) { + return a.type > b.type; + } - for (size_t i=0; i SeamPlacer::overhang_distance_tolerance_factor * a.perimeter.flow_width || + b.overhang > SeamPlacer::overhang_distance_tolerance_factor * b.perimeter.flow_width) { + return a.overhang < b.overhang; + } + + // prefer hidden points (more than 1 mm inside) + if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) { + return true; + } + if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) { + return false; + } + + if (setup == SeamPosition::spRear) { + if (a.position.y() == b.position.y()) { + return a.position.x() > b.position.x(); + } + 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.visibility + + SeamPlacer::angle_importance * compute_angle_penalty(a.local_ccw_angle) + + distance_penalty_a; + float penalty_b = b.visibility + + SeamPlacer::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 > SeamPlacer::overhang_distance_tolerance_factor * a.perimeter.flow_width || + b.overhang > SeamPlacer::overhang_distance_tolerance_factor * b.perimeter.flow_width) { + return a.overhang < b.overhang; + } + + // prefer hidden points (more than 1 mm inside) + if (a.embedded_distance < -1.0f && b.embedded_distance > -1.0f) { + return true; + } + if (b.embedded_distance < -1.0f && a.embedded_distance > -1.0f) { + return false; + } + + if (setup == SeamPosition::spRandom) { + return true; + } + + if (setup == SeamPosition::spRear) { + if (a.position.y() == b.position.y()) { + return a.position.x() > b.position.x(); + } + return a.position.y() > b.position.y(); + } + + float penalty_a = a.visibility + + SeamPlacer::angle_importance * compute_angle_penalty(a.local_ccw_angle); + float penalty_b = b.visibility + + SeamPlacer::angle_importance * compute_angle_penalty(b.local_ccw_angle); + + return penalty_a <= penalty_b || penalty_a - penalty_b < SeamPlacer::seam_align_score_tolerance; + } + + bool are_similar(const SeamCandidate &a, const SeamCandidate &b) const { + return is_first_not_much_worse(a, b) && is_first_not_much_worse(b, a); + } + + float compute_angle_penalty(float ccw_angle) const { + // This function is used: + // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x))) + // looks scary, but it is gaussian combined with sigmoid, + // so that concave points have much smaller penalty over convex ones + // https://github.com/prusa3d/PrusaSlicer/tree/master/doc/seam_placement/corner_penalty_function.png + return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) + + 1.0f / (2 + std::exp(-ccw_angle)); + } +}; + +#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, -comparator.compute_angle_penalty(point.local_ccw_angle)); + max_weight = std::max(max_weight, -comparator.compute_angle_penalty(point.local_ccw_angle)); + + } + + std::string visiblity_file_name = debug_out_path( + ("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, + -comparator.compute_angle_penalty(point.local_ccw_angle)); + std::string weight_fill = "rgb(" + std::to_string(weight_color.x()) + "," + std::to_string(weight_color.y()) + + "," + + std::to_string(weight_color.z()) + ")"; + weight_svg.draw(scaled(Vec2f(point.position.head<2>())), weight_fill); + + Vec3i overhang_color = value_to_rgbi(-0.5, 0.5, std::clamp(point.overhang, -0.5f, 0.5f)); + std::string overhang_fill = "rgb(" + std::to_string(overhang_color.x()) + "," + + std::to_string(overhang_color.y()) + + "," + + std::to_string(overhang_color.z()) + ")"; + overhangs_svg.draw(scaled(Vec2f(point.position.head<2>())), overhang_fill); + } + } +} +#endif + +// Pick best seam point based on the given comparator +void pick_seam_point(std::vector &perimeter_points, size_t start_index, + const SeamComparator &comparator) { + size_t end_index = perimeter_points[start_index].perimeter.end_index; + + size_t seam_index = start_index; + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index])) { + seam_index = index; + } + } + perimeter_points[start_index].perimeter.seam_index = seam_index; +} + +size_t pick_nearest_seam_point_index(const std::vector &perimeter_points, size_t start_index, + const Vec2f &preffered_location) { + size_t end_index = perimeter_points[start_index].perimeter.end_index; + SeamComparator comparator { spNearest }; + + size_t seam_index = start_index; + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.is_first_better(perimeter_points[index], perimeter_points[seam_index], preffered_location)) { + seam_index = index; + } + } + return seam_index; +} + +// picks random seam point uniformly, respecting enforcers blockers and overhang avoidance. +void pick_random_seam_point(const std::vector &perimeter_points, size_t start_index) { + SeamComparator comparator { spRandom }; + + // algorithm keeps a list of viable points and their lengths. If it finds a point + // that is much better than the viable_example_index (e.g. better type, no overhang; see is_first_not_much_worse) + // then it throws away stored lists and starts from start + // in the end, the list should contain points with same type (Enforced > Neutral > Blocked) and also only those which are not + // big overhang. + size_t viable_example_index = start_index; + size_t end_index = perimeter_points[start_index].perimeter.end_index; + struct Viable { + // Candidate seam point index. + size_t index; + float edge_length; + Vec3f edge; + }; + std::vector viables; + + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.are_similar(perimeter_points[index], perimeter_points[viable_example_index])) { + // index ok, push info into viables + Vec3f edge_to_next { perimeter_points[index == end_index ? start_index : index + 1].position + - perimeter_points[index].position }; + float dist_to_next = edge_to_next.norm(); + viables.push_back( { index, dist_to_next, edge_to_next }); + } else if (comparator.is_first_not_much_worse(perimeter_points[viable_example_index], + perimeter_points[index])) { + // index is worse then viable_example_index, skip this point + } else { + // index is better than viable example index, update example, clear gathered info, start again + // clear up all gathered info, start from scratch, update example index + viable_example_index = index; + viables.clear(); + + Vec3f edge_to_next = (perimeter_points[index == end_index ? start_index : index + 1].position + - perimeter_points[index].position); + float dist_to_next = edge_to_next.norm(); + viables.push_back( { index, dist_to_next, edge_to_next }); + } + } + + // now pick random point from the stored options + float len_sum = std::accumulate(viables.begin(), viables.end(), 0.0f, [](const float acc, const Viable &v) { + return acc + v.edge_length; + }); + float picked_len = len_sum * (rand() / (float(RAND_MAX) + 1)); + + size_t point_idx = 0; + while (picked_len - viables[point_idx].edge_length > 0) { + picked_len = picked_len - viables[point_idx].edge_length; + point_idx++; + } + + Perimeter &perimeter = perimeter_points[start_index].perimeter; + perimeter.seam_index = viables[point_idx].index; + perimeter.final_seam_position = perimeter_points[perimeter.seam_index].position + + viables[point_idx].edge.normalized() * picked_len; + perimeter.finalized = true; +} + +class PerimeterDistancer { + std::vector lines; + AABBTreeIndirect::Tree<2, double> tree; + +public: + PerimeterDistancer(const Layer *layer) { + static const float eps = float(scale_(layer->object()->config().slice_closing_radius.value)); + // merge with offset + ExPolygons merged = layer->merged(eps); + // ofsset back + ExPolygons layer_outline = offset_ex(merged, -eps); + 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; + 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) { + bool layer_has_multiple_loops = + layers[layer_idx].points[0].perimeter.end_index + < layers[layer_idx].points.size() - 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); + } + + if (layer_has_multiple_loops) { // 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); + } + } + + 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 std::pair &prev_point_index, + const size_t layer_idx, const float slice_z, + const SeamPlacerImpl::SeamComparator &comparator) const { + using namespace SeamPlacerImpl; + + const SeamCandidate &last_point = layers[prev_point_index.first].points[prev_point_index.second]; + + Vec3f projected_position { last_point.position.x(), last_point.position.y(), slice_z }; + std::vector nearby_points_indices = find_nearby_points(*layers[layer_idx].points_tree, projected_position, + SeamPlacer::seam_align_tolerable_dist); + + if (nearby_points_indices.empty()) { + return {}; + } + + 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 * SeamPlacer::seam_align_tolerable_dist)) { + 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, + std::optional> &out_best_moved_seam, size_t &out_moved_seams_count) const { + out_best_moved_seam.reset(); + out_moved_seams_count = 0; + const std::vector &layers = m_seam_per_object.find(po)->second.layers; + int layer_idx = start_seam.first; + int seam_index = start_seam.second; + + //initialize searching for seam string - cluster of nearby seams on previous and next layers + int skips = SeamPlacer::seam_align_tolerable_skips / 2; + int next_layer = layer_idx + 1; + std::pair prev_point_index = start_seam; + std::vector> seam_string { start_seam }; + + //find seams or potential seams in forward direction; there is a budget of skips allowed + while (skips >= 0 && next_layer < int(layers.size())) { + auto maybe_next_seam = find_next_seam_in_layer(layers, prev_point_index, next_layer, + float(po->get_layer(next_layer)->slice_z), comparator); + if (maybe_next_seam.has_value()) { + + std::pair next_seam_coords = maybe_next_seam.value(); + const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second]; + bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second; + out_moved_seams_count += is_moved; + if (is_moved && (!out_best_moved_seam.has_value() || + comparator.is_first_better(next_seam, + layers[out_best_moved_seam.value().first].points[out_best_moved_seam.value().second]))) { + out_best_moved_seam = { next_seam_coords }; + } + + seam_string.push_back(maybe_next_seam.value()); + prev_point_index = seam_string.back(); + //String added, prev_point_index updated + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer++; + } + + //do additional check in back direction + next_layer = layer_idx - 1; + skips = SeamPlacer::seam_align_tolerable_skips / 2; + prev_point_index = std::pair(layer_idx, seam_index); + while (skips >= 0 && next_layer >= 0) { + auto maybe_next_seam = find_next_seam_in_layer(layers, prev_point_index, next_layer, + float(po->get_layer(next_layer)->slice_z), comparator); + if (maybe_next_seam.has_value()) { + + std::pair next_seam_coords = maybe_next_seam.value(); + const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second]; + bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second; + out_moved_seams_count += is_moved; + if (is_moved && (!out_best_moved_seam.has_value() || + comparator.is_first_better(next_seam, + layers[out_best_moved_seam.value().first].points[out_best_moved_seam.value().second]))) { + out_best_moved_seam = { next_seam_coords }; + } + + seam_string.push_back(maybe_next_seam.value()); + prev_point_index = seam_string.back(); + //String added, prev_point_index updated + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer--; + } + + return seam_string; +} + +// clusters already chosen seam points into strings across multiple layers, and then +// aligns the strings via polynomial fit +// Does not change the positions of the SeamCandidates themselves, instead stores +// the new aligned position into the shared Perimeter structure of each perimeter +// Note that this position does not necesarilly lay on the perimeter. +void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) { + using namespace SeamPlacerImpl; + + // Prepares Debug files for writing. +#ifdef DEBUG_FILES + Slic3r::CNumericLocalesSetter locales_setter; + auto clusters_f = debug_out_path("seam_clusters.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; + } + auto aligned_f = debug_out_path("aligned_clusters.obj"); + FILE *aligns = boost::nowide::fopen(aligned_f.c_str(), "w"); + if (aligns == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << clusters_f << " for writing"; + return; + } +#endif + + //gather vector of all seams on the print_object - pair of layer_index and seam__index within that layer + const std::vector &layers = m_seam_per_object[po].layers; + std::vector> seams; + for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) { + const std::vector &layer_perimeter_points = layers[layer_idx].points; + size_t current_point_index = 0; + while (current_point_index < layer_perimeter_points.size()) { + seams.emplace_back(layer_idx, layer_perimeter_points[current_point_index].perimeter.seam_index); + current_point_index = layer_perimeter_points[current_point_index].perimeter.end_index + 1; + } + } + + //sort them before alignment. Alignment is sensitive to initializaion, this gives it better chance to choose something nice + std::sort(seams.begin(), seams.end(), + [&comparator, &layers](const std::pair &left, + const std::pair &right) { + return comparator.is_first_better(layers[left.first].points[left.second], + layers[right.first].points[right.second]); + } + ); + + //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment + // Keeping the vectors outside, so with a bit of luck they will not get reallocated after couple of for loop iterations. + std::vector> seam_string; + std::vector> alternative_seam_string; + std::vector observations; + std::vector observation_points; + std::vector weights; + + 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 { + std::optional> best_moved_seam; + size_t moved_seams_count; + seam_string = this->find_seam_string(po, { layer_idx, seam_index }, comparator, best_moved_seam, + moved_seams_count); + if (best_moved_seam.has_value()) { + size_t alternative_moved_seams_count; + alternative_seam_string = this->find_seam_string(po, best_moved_seam.value(), comparator, + best_moved_seam, alternative_moved_seams_count); + if (alternative_seam_string.size() >= SeamPlacer::seam_align_minimum_string_seams && + alternative_moved_seams_count < moved_seams_count) { + seam_string = std::move(alternative_seam_string); + // finish loop. but repeat the alignment for the current seam, since it could be skipped due to alternative path being aligned. + global_index--; + } + } + + if (seam_string.size() < seam_align_minimum_string_seams) { + //string NOT long enough to be worth aligning, skip + continue; + } + + // String is long enough, all string seams and potential string seams gathered, now do the alignment + //sort by layer index + std::sort(seam_string.begin(), seam_string.end(), + [](const std::pair &left, const std::pair &right) { + return left.first < right.first; + }); + + // gather all positions of seams and their weights (weights are derived as negative penalty, they are made positive in next step) + observations.resize(seam_string.size()); + observation_points.resize(seam_string.size()); + weights.resize(seam_string.size()); + + //gather points positions and weights + // The algorithm uses only angle to compute penalty, to enforce snapping to sharp corners, if they are present + // after several experiments approach that gives best results is to snap the weight to one for sharp corners, and + // leave it small for others. However, this can result in non-smooth line over area with a lot of unaligned sharp corners. + for (size_t index = 0; index < seam_string.size(); ++index) { + Vec3f pos = layers[seam_string[index].first].points[seam_string[index].second].position; + observations[index] = pos.head<2>(); + observation_points[index] = pos.z(); + weights[index] = + (comparator.compute_angle_penalty( + layers[seam_string[index].first].points[seam_string[index].second].local_ccw_angle) + < comparator.compute_angle_penalty(0.4f * float(PI))) ? 1.0f : 0.1f; + } + + // Curve Fitting + size_t number_of_segments = std::max(size_t(1), + size_t(observations.size() / SeamPlacer::seam_align_seams_per_segment)); + auto curve = Geometry::fit_cubic_bspline(observations, observation_points, weights, number_of_segments); + + // Do alignment - compute fitted point for each point in the string from its Z coord, and store the position into + // Perimeter structure of the point; also set flag aligned to true + for (size_t index = 0; index < seam_string.size(); ++index) { + const auto &pair = seam_string[index]; + const float t = weights[index]; + Vec3f current_pos = layers[pair.first].points[pair.second].position; + Vec2f fitted_pos = curve.get_fitted_value(current_pos.z()); + + //interpolate between current and fitted position, prefer current pos for large weights. + Vec3f final_position = t * current_pos + (1 - t) * to_3d(fitted_pos, current_pos.z()); + + Perimeter &perimeter = layers[pair.first].points[pair.second].perimeter; + perimeter.seam_index = pair.second; + perimeter.final_seam_position = final_position; + perimeter.finalized = true; + } + +#ifdef DEBUG_FILES + auto randf = []() { + return float(rand()) / float(RAND_MAX); + }; + Vec3f color { randf(), randf(), randf() }; + for (size_t i = 0; i < seam_string.size(); ++i) { + auto orig_seam = layers[seam_string[i].first].points[seam_string[i].second]; + fprintf(clusters, "v %f %f %f %f %f %f \n", orig_seam.position[0], + orig_seam.position[1], + orig_seam.position[2], color[0], color[1], + color[2]); + } + + color = Vec3f { randf(), randf(), randf() }; + for (size_t i = 0; i < seam_string.size(); ++i) { + const Perimeter &perimeter = layers[seam_string[i].first].points[seam_string[i].second].perimeter; + fprintf(aligns, "v %f %f %f %f %f %f \n", perimeter.final_seam_position[0], + perimeter.final_seam_position[1], + perimeter.final_seam_position[2], color[0], color[1], + color[2]); + } +#endif + } + } + +#ifdef DEBUG_FILES + fclose(clusters); + fclose(aligns); #endif } +void SeamPlacer::init(const Print &print) { + using namespace SeamPlacerImpl; + m_seam_per_object.clear(); + for (const PrintObject *po : print.objects()) { -std::optional SeamHistory::get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb) -{ - assert(layer_id >= m_layer_id || layer_id == 0); - if (layer_id != m_layer_id) { - // Get seam was called for different layer than last time. - if (layer_id == 0) // seq printing - m_data_this_layer.clear(); - m_data_last_layer = m_data_this_layer; - m_data_this_layer.clear(); - m_layer_id = layer_id; - } + SeamPosition configured_seam_preference = po->config().seam_position.value; + SeamComparator comparator { configured_seam_preference }; - std::optional out; + GlobalModelInfo global_model_info { }; + gather_enforcers_blockers(global_model_info, po); - auto seams_it = m_data_last_layer.find(po); - if (seams_it == m_data_last_layer.end()) - return out; - - const std::vector& seam_data_po = seams_it->second; - - // Find a bounding-box on the last layer that is close to one we see now. - double min_score = std::numeric_limits::max(); - for (const SeamPoint& sp : seam_data_po) { - const BoundingBox& bb = sp.m_island_bb; - - if (! bb.overlap(island_bb)) { - // This bb does not even overlap. It is likely unrelated. - continue; + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { + compute_global_occlusion(global_model_info, po); } - double score = std::pow(bb.min(0) - island_bb.min(0), 2.) - + std::pow(bb.min(1) - island_bb.min(1), 2.) - + std::pow(bb.max(0) - island_bb.max(0), 2.) - + std::pow(bb.max(1) - island_bb.max(1), 2.); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather_seam_candidates: start"; + gather_seam_candidates(po, global_model_info, configured_seam_preference); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: gather_seam_candidates: end"; - if (score < min_score) { - min_score = score; - out = sp.m_pos; + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_candidates_visibility : start"; + calculate_candidates_visibility(po, global_model_info); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_candidates_visibility : end"; } + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_overhangs and layer embdedding : start"; + calculate_overhangs_and_layer_embedding(po); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_overhangs and layer embdedding: end"; + + if (configured_seam_preference != spNearest) { // For spNearest, the seam is picked in the place_seam method with actual nozzle position information + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: pick_seam_point : start"; + //pick seam point + std::vector &layers = m_seam_per_object[po].layers; + tbb::parallel_for(tbb::blocked_range(0, layers.size()), + [&layers, configured_seam_preference, comparator](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + std::vector &layer_perimeter_points = layers[layer_idx].points; + for (size_t current = 0; current < layer_perimeter_points.size(); + current = layer_perimeter_points[current].perimeter.end_index + 1) + if (configured_seam_preference == spRandom) + pick_random_seam_point(layer_perimeter_points, current); + else + pick_seam_point(layer_perimeter_points, current, comparator); + } + }); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: pick_seam_point : end"; + } + + if (configured_seam_preference == spAligned) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: align_seam_points : start"; + align_seam_points(po, comparator); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: align_seam_points : end"; + } + +#ifdef DEBUG_FILES + debug_export_points(m_seam_per_object[po].layers, po->bounding_box(), comparator); +#endif + } +} + +void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, + const Point &last_pos) const { + using namespace SeamPlacerImpl; + const PrintObject *po = layer->object(); + // Must not be called with supprot layer. + assert(dynamic_cast(layer) == nullptr); + // Object layer IDs are incremented by the number of raft layers. + assert(layer->id() >= po->slicing_parameters().raft_layers()); + const size_t layer_index = layer->id() - po->slicing_parameters().raft_layers(); + const double unscaled_z = layer->slice_z; + + const PrintObjectSeamData::LayerSeams &layer_perimeters = + m_seam_per_object.find(layer->object())->second.layers[layer_index]; + + // Find the closest perimeter in the SeamPlacer to the first point of this loop. + size_t closest_perimeter_point_index; + { + const Point &fp = loop.first_point(); + Vec2f unscaled_p = unscaled(fp); + closest_perimeter_point_index = find_closest_point(*layer_perimeters.points_tree.get(), + to_3d(unscaled_p, float(unscaled_z))); } - return out; -} + Vec3f seam_position; + size_t seam_index; + if (const Perimeter &perimeter = layer_perimeters.points[closest_perimeter_point_index].perimeter ; + perimeter.finalized) { + seam_position = perimeter.final_seam_position; + seam_index = perimeter.seam_index; + } else { + seam_index = + po->config().seam_position == spNearest ? + pick_nearest_seam_point_index(layer_perimeters.points, perimeter.start_index, + unscaled(last_pos)) : + perimeter.seam_index; + seam_position = layer_perimeters.points[seam_index].position; + } + Point seam_point = Point::new_scale(seam_position.x(), seam_position.y()); + if (const SeamCandidate &perimeter_point = layer_perimeters.points[seam_index]; + (po->config().seam_position == spNearest || po->config().seam_position == spAligned) && + loop.role() == ExtrusionRole::erPerimeter && //Hopefully internal perimeter + (seam_position - perimeter_point.position).squaredNorm() < 4.0f && // seam is on perimeter point + perimeter_point.local_ccw_angle < -EPSILON // In concave angles + ) { // In this case, we are at internal perimeter, where the external perimeter has seam in concave angle. We want to align + // the internal seam into the concave corner, and not on the perpendicular projection on the closest edge (which is what the split_at function does) + size_t index_of_prev = + seam_index == perimeter_point.perimeter.start_index ? + perimeter_point.perimeter.end_index : + seam_index - 1; + size_t index_of_next = + seam_index == perimeter_point.perimeter.end_index ? + perimeter_point.perimeter.start_index : + seam_index + 1; -void SeamHistory::add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb) -{ - m_data_this_layer[po].push_back({pos, island_bb});; -} + Vec2f dir_to_middle = + ((perimeter_point.position - layer_perimeters.points[index_of_prev].position).head<2>().normalized() + + (perimeter_point.position - layer_perimeters.points[index_of_next].position).head<2>().normalized()) + * 0.5; + auto [_, projected_point] = loop.get_closest_path_and_point(seam_point, true); + //get closest projected point, determine depth of the seam point. + float depth = (float) unscale(Point(seam_point - projected_point)).norm(); + float angle_factor = cos(-perimeter_point.local_ccw_angle / 2.0f); // There are some nice geometric identities in determination of the correct depth of new seam point. + //overshoot the target depth, in concave angles it will correctly snap to the corner; TODO: find out why such big overshoot is needed. + Vec2f final_pos = perimeter_point.position.head<2>() + (1.4142 * depth / angle_factor) * dir_to_middle; + seam_point = Point::new_scale(final_pos.x(), final_pos.y()); + } - -void SeamHistory::clear() -{ - m_layer_id = 0; - m_data_last_layer.clear(); - m_data_this_layer.clear(); -} - + if (!loop.split_at_vertex(seam_point)) { + // The point is not in the original loop. + // Insert it. + loop.split_at(seam_point, true); + } } + +} // namespace Slic3r diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 0d72939b2..7dbe7bcb3 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,173 @@ class PrintObject; class ExtrusionLoop; class Print; class Layer; -namespace EdgeGrid { class Grid; } + +namespace EdgeGrid { +class Grid; +} + +namespace SeamPlacerImpl { -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(); +// ************ 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)); +} +// *************************** -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; +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; +}; + +//Struct over which all processing of perimeters is done. For each perimeter point, its respective candidate is created, +// then all the needed attributes are computed and finally, for each perimeter one point is chosen as seam. +// This seam position can be then further aligned +struct SeamCandidate { + SeamCandidate(const Vec3f &pos, Perimeter &perimeter, + float local_ccw_angle, + EnforcedBlockedSeamPoint type) : + position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), embedded_distance(0.0f), local_ccw_angle( + local_ccw_angle), type(type), central_enforcer(false) { + } + const Vec3f position; + // pointer to Perimeter loop of this point. It is shared across all points of the loop + Perimeter &perimeter; + float visibility; + float overhang; + // distance inside the merged layer regions, for detecting perimeter points which are hidden indside the print (e.g. multimaterial join) + // Negative sign means inside the print, comes from EdgeGrid structure + float embedded_distance; + float local_ccw_angle; + EnforcedBlockedSeamPoint type; + bool central_enforcer; //marks this candidate as central point of enforced segment on the perimeter - important for alignment +}; + +struct FaceVisibilityInfo { + float visibility; +}; + +struct SeamCandidateCoordinateFunctor { + SeamCandidateCoordinateFunctor(const std::vector &seam_candidates) : + seam_candidates(seam_candidates) { + } + const std::vector &seam_candidates; + float operator()(size_t index, size_t dim) const { + return seam_candidates[index].position[dim]; + } +}; +} // namespace SeamPlacerImpl + +struct PrintObjectSeamData +{ + using SeamCandidatesTree = KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>; + + struct LayerSeams + { + Slic3r::deque perimeters; + std::vector points; + std::unique_ptr points_tree; + }; + // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter + std::vector layers; + // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD + // tree of all points of the given layer + + void clear() + { + layers.clear(); + } +}; class SeamPlacer { public: - void init(const Print& print); + static constexpr size_t raycasting_decimation_target_triangle_count = 10000; + // for subdivision, both following criteria are considered, and the one with less resulting triangles is used + static constexpr size_t raycasting_subdivision_target_triangle_count = 20000; + static constexpr float raycasting_subdivision_target_length = 2.0f; + //square of number of rays per triangle + static constexpr size_t sqr_rays_per_triangle = 7; - // When perimeters are printed, first call this function with the respective - // external perimeter. SeamPlacer will find a location for its seam and remember it. - // Subsequent calls to get_seam will return this position. + // arm length used during angles computation + static constexpr float polygon_local_angles_arm_distance = 0.1f; - 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); + // max tolerable distance from the previous layer is overhang_distance_tolerance_factor * flow_width + static constexpr float overhang_distance_tolerance_factor = 0.5f; - 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; + // determines angle importance compared to visibility ( neutral value is 1.0f. ) + static constexpr float angle_importance = 0.7f; + + // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer + static constexpr float enforcer_blocker_distance_tolerance = 0.35f; + // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size + static constexpr float enforcer_oversampling_distance = 0.2f; + + // 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.25f; + // seam_align_tolerable_dist - if next layer closes point is too far away, break string + static constexpr float seam_align_tolerable_dist = 1.0f; + // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. + // this param limits the number of allowed skips + static constexpr size_t seam_align_tolerable_skips = 4; + // minimum number of seams needed in cluster to make alignment happen + static constexpr size_t seam_align_minimum_string_seams = 6; + // points covered by spline; determines number of splines for the given string + static constexpr size_t seam_align_seams_per_segment = 8; + + //The following data structures hold all perimeter points for all PrintObject. + std::unordered_map m_seam_per_object; + + void init(const Print &print); + + void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const; private: - - // When given an external perimeter (!), returns the seam. - Point calculate_seam(const Layer& layer, const SeamPosition seam_position, - const ExtrusionLoop& loop, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid, Point last_pos, bool prefer_nearest); - - struct CustomTrianglesPerLayer { - Polygons polys; - TreeType tree; - }; - - // Just a cache to save some lookups. - const Layer* m_last_layer_po = nullptr; - coordf_t m_last_print_z = -1.; - const PrintObject* m_last_po = nullptr; - - struct SeamPoint { - Point pt; - bool precalculated = false; - bool external = false; - const Layer* layer = nullptr; - SeamPosition seam_position; - const PrintObject* po = nullptr; - }; - std::vector m_plan; - size_t m_plan_idx; - - std::vector> m_enforcers; - std::vector> m_blockers; - std::vector m_po_list; - - //std::map m_last_seam_position; - SeamHistory m_seam_history; - - // Get indices of points inside enforcers and blockers. - void get_enforcers_and_blockers(size_t layer_id, - const Polygon& polygon, - size_t po_id, - std::vector& enforcers_idxs, - std::vector& blockers_idxs) const; - - // Apply penalties to points inside enforcers/blockers. - void apply_custom_seam(const Polygon& polygon, size_t po_id, - std::vector& penalties, - const std::vector& lengths, - int layer_id, SeamPosition seam_position) const; - - // Return random point of a polygon. The distribution will be uniform - // along the contour and account for enforcers and blockers. - Point get_random_seam(size_t layer_idx, const Polygon& polygon, size_t po_id, - bool* saw_custom = nullptr) const; - - // Is there any enforcer/blocker on this layer? - bool is_custom_seam_on_layer(size_t layer_id, size_t po_idx) const { - return is_custom_enforcer_on_layer(layer_id, po_idx) - || is_custom_blocker_on_layer(layer_id, po_idx); - } - - bool is_custom_enforcer_on_layer(size_t layer_id, size_t po_idx) const { - return (! m_enforcers.at(po_idx).empty() && ! m_enforcers.at(po_idx)[layer_id].polys.empty()); - } - - bool is_custom_blocker_on_layer(size_t layer_id, size_t po_idx) const { - return (! m_blockers.at(po_idx).empty() && ! m_blockers.at(po_idx)[layer_id].polys.empty()); - } + void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info, + const SeamPosition configured_seam_preference); + void calculate_candidates_visibility(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info); + void calculate_overhangs_and_layer_embedding(const PrintObject *po); + void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator); + std::vector> find_seam_string(const PrintObject *po, + std::pair start_seam, + const SeamPlacerImpl::SeamComparator &comparator, + std::optional> &out_best_moved_seam, + size_t& out_moved_seams_count) const; + std::optional> find_next_seam_in_layer( + const std::vector &layers, + const std::pair &prev_point_index, + const size_t layer_idx, const float slice_z, + 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 000000000..98c6f8bb2 --- /dev/null +++ b/src/libslic3r/Geometry/Bicubic.hpp @@ -0,0 +1,291 @@ +#ifndef BICUBIC_HPP +#define BICUBIC_HPP + +#include +#include +#include + +#include + +namespace Slic3r { + +namespace Geometry { + +namespace BicubicInternal { +// Linear kernel, to be able to test cubic methods with hat kernels. +template +struct LinearKernel +{ + typedef T FloatType; + + static T a00() { + return T(0.); + } + static T a01() { + return T(0.); + } + static T a02() { + return T(0.); + } + static T a03() { + return T(0.); + } + static T a10() { + return T(1.); + } + static T a11() { + return T(-1.); + } + static T a12() { + return T(0.); + } + static T a13() { + return T(0.); + } + static T a20() { + return T(0.); + } + static T a21() { + return T(1.); + } + static T a22() { + return T(0.); + } + static T a23() { + return T(0.); + } + static T a30() { + return T(0.); + } + static T a31() { + return T(0.); + } + static T a32() { + return T(0.); + } + static T a33() { + return T(0.); + } +}; + +// Interpolation kernel aka Catmul-Rom aka Keyes kernel. +template +struct CubicCatmulRomKernel +{ + typedef T FloatType; + + static T a00() { + return 0; + } + static T a01() { + return T( -0.5); + } + static T a02() { + return T( 1.); + } + static T a03() { + return T( -0.5); + } + static T a10() { + return T( 1.); + } + static T a11() { + return 0; + } + static T a12() { + return T( -5. / 2.); + } + static T a13() { + return T( 3. / 2.); + } + static T a20() { + return 0; + } + static T a21() { + return T( 0.5); + } + static T a22() { + return T( 2.); + } + static T a23() { + return T( -3. / 2.); + } + static T a30() { + return 0; + } + static T a31() { + return 0; + } + static T a32() { + return T( -0.5); + } + static T a33() { + return T( 0.5); + } +}; + +// B-spline kernel +template +struct CubicBSplineKernel +{ + typedef T FloatType; + + static T a00() { + return T( 1. / 6.); + } + static T a01() { + return T( -3. / 6.); + } + static T a02() { + return T( 3. / 6.); + } + static T a03() { + return T( -1. / 6.); + } + static T a10() { + return T( 4. / 6.); + } + static T a11() { + return 0; + } + static T a12() { + return T( -6. / 6.); + } + static T a13() { + return T( 3. / 6.); + } + static T a20() { + return T( 1. / 6.); + } + static T a21() { + return T( 3. / 6.); + } + static T a22() { + return T( 3. / 6.); + } + static T a23() { + return T( -3. / 6.); + } + static T a30() { + return 0; + } + static T a31() { + return 0; + } + static T a32() { + return 0; + } + static T a33() { + return T( 1. / 6.); + } +}; + +template +inline T clamp(T a, T lower, T upper) + { + return (a < lower) ? lower : + (a > upper) ? upper : a; +} +} + +template +struct CubicKernelWrapper +{ + typedef typename Kernel::FloatType FloatType; + + static constexpr size_t kernel_span = 4; + + static FloatType kernel(FloatType x) + { + x = fabs(x); + if (x >= (FloatType) 2.) + return 0.0f; + if (x <= (FloatType) 1.) { + FloatType x2 = x * x; + FloatType x3 = x2 * x; + return Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3; + } + assert(x > (FloatType )1. && x < (FloatType )2.); + x -= (FloatType) 1.; + FloatType x2 = x * x; + FloatType x3 = x2 * x; + return Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3; + } + + static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x) + { + const FloatType x2 = x * x; + const FloatType x3 = x * x * x; + return f0 * (Kernel::a00() + Kernel::a01() * x + Kernel::a02() * x2 + Kernel::a03() * x3) + + f1 * (Kernel::a10() + Kernel::a11() * x + Kernel::a12() * x2 + Kernel::a13() * x3) + + f2 * (Kernel::a20() + Kernel::a21() * x + Kernel::a22() * x2 + Kernel::a23() * x3) + + f3 * (Kernel::a30() + Kernel::a31() * x + Kernel::a32() * x2 + Kernel::a33() * x3); + } +}; + +// Linear splines +template +using LinearKernel = CubicKernelWrapper>; + +// Catmul-Rom splines +template +using CubicCatmulRomKernel = CubicKernelWrapper>; + +// Cubic B-splines +template +using CubicBSplineKernel = CubicKernelWrapper>; + +template +static typename KernelWrapper::FloatType cubic_interpolate(const Eigen::ArrayBase &F, + const typename KernelWrapper::FloatType pt) { + typedef typename KernelWrapper::FloatType T; + const int w = int(F.size()); + const int ix = (int) floor(pt); + const T s = pt - T( ix); + + if (ix > 1 && ix + 2 < w) { + // Inside the fully interpolated region. + return KernelWrapper::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s); + } + // Transition region. Extend with a constant function. + auto f = [&F, w](T x) { + return F[BicubicInternal::clamp(x, 0, w - 1)]; + }; + return KernelWrapper::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s); +} + +template +static float bicubic_interpolate(const Eigen::MatrixBase &F, + const Eigen::Matrix &pt) { + typedef typename Kernel::FloatType T; + const int w = F.cols(); + const int h = F.rows(); + const int ix = (int) floor(pt[0]); + const int iy = (int) floor(pt[1]); + const T s = pt[0] - T( ix); + const T t = pt[1] - T( iy); + + if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) { + // Inside the fully interpolated region. + return Kernel::interpolate( + Kernel::interpolate(F(ix - 1, iy - 1), F(ix, iy - 1), F(ix + 1, iy - 1), F(ix + 2, iy - 1), s), + Kernel::interpolate(F(ix - 1, iy), F(ix, iy), F(ix + 1, iy), F(ix + 2, iy), s), + Kernel::interpolate(F(ix - 1, iy + 1), F(ix, iy + 1), F(ix + 1, iy + 1), F(ix + 2, iy + 1), s), + Kernel::interpolate(F(ix - 1, iy + 2), F(ix, iy + 2), F(ix + 1, iy + 2), F(ix + 2, iy + 2), s), t); + } + // Transition region. Extend with a constant function. + auto f = [&F, w, h](int x, int y) { + return F(BicubicInternal::clamp(x, 0, w - 1), BicubicInternal::clamp(y, 0, h - 1)); + }; + return Kernel::interpolate( + Kernel::interpolate(f(ix - 1, iy - 1), f(ix, iy - 1), f(ix + 1, iy - 1), f(ix + 2, iy - 1), s), + Kernel::interpolate(f(ix - 1, iy), f(ix, iy), f(ix + 1, iy), f(ix + 2, iy), s), + Kernel::interpolate(f(ix - 1, iy + 1), f(ix, iy + 1), f(ix + 1, iy + 1), f(ix + 2, iy + 1), s), + Kernel::interpolate(f(ix - 1, iy + 2), f(ix, iy + 2), f(ix + 1, iy + 2), f(ix + 2, iy + 2), s), t); +} + +} //namespace Geometry + +} // namespace Slic3r + +#endif /* BICUBIC_HPP */ diff --git a/src/libslic3r/Geometry/Curves.hpp b/src/libslic3r/Geometry/Curves.hpp new file mode 100644 index 000000000..6542cb706 --- /dev/null +++ b/src/libslic3r/Geometry/Curves.hpp @@ -0,0 +1,205 @@ +#ifndef SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ +#define SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ + +#include "libslic3r/Point.hpp" +#include "Bicubic.hpp" + +#include + +//#define LSQR_DEBUG + +namespace Slic3r { +namespace Geometry { + +template +struct PolynomialCurve { + Eigen::MatrixXf coefficients; + + Vec3f get_fitted_value(const NumberType value) const { + auto result = Vec::Zero(); + size_t order = this->coefficients.rows() - 1; + auto x = NumberType(1.); + for (size_t index = 0; index < order + 1; ++index, x *= value) + result += x * this->coefficients.col(index); + return result; + } +}; + +//https://towardsdatascience.com/least-square-polynomial-CURVES-using-c-eigen-package-c0673728bd01 +template +PolynomialCurve fit_polynomial(const std::vector> &observations, + const std::vector &observation_points, + const std::vector &weights, size_t order) { + // check to make sure inputs are correct + size_t cols = order + 1; + assert(observation_points.size() >= cols); + assert(observation_points.size() == weights.size()); + assert(observations.size() == weights.size()); + + Eigen::MatrixXf data_points(Dimension, observations.size()); + Eigen::MatrixXf T(observations.size(), cols); + for (size_t i = 0; i < weights.size(); ++i) { + auto squared_weight = sqrt(weights[i]); + data_points.col(i) = observations[i] * squared_weight; + // Populate the matrix + auto x = squared_weight; + auto c = observation_points[i]; + for (size_t j = 0; j < cols; ++j, x *= c) + T(i, j) = x; + } + + const auto QR = T.householderQr(); + Eigen::MatrixXf coefficients(Dimension, cols); + // Solve for linear least square fit + for (size_t dim = 0; dim < Dimension; ++dim) { + coefficients.row(dim) = QR.solve(data_points.row(dim).transpose()); + } + + return {std::move(coefficients)}; +} + +template +struct PiecewiseFittedCurve { + using Kernel = KernelType; + + Eigen::MatrixXf coefficients; + NumberType start; + NumberType segment_size; + size_t endpoints_level_of_freedom; + + Vec get_fitted_value(const NumberType &observation_point) const { + Vec result = Vec::Zero(); + + //find corresponding segment index; expects kernels to be centered + int middle_right_segment_index = floor((observation_point - start) / segment_size); + //find index of first segment that is affected by the point i; this can be deduced from kernel_span + int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1; + for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span); + segment_index++) { + NumberType segment_start = start + segment_index * segment_size; + NumberType normalized_segment_distance = (segment_start - observation_point) / segment_size; + + int parameter_index = segment_index + endpoints_level_of_freedom; + parameter_index = std::clamp(parameter_index, 0, int(coefficients.cols()) - 1); + result += Kernel::kernel(normalized_segment_distance) * coefficients.col(parameter_index); + } + return result; + } +}; + +// observations: data to be fitted by the curve +// observation points: growing sequence of points where the observations were made. +// In other words, for function f(x) = y, observations are y0...yn, and observation points are x0...xn +// weights: how important the observation is +// segments_count: number of segments inside the valid length of the curve +// endpoints_level_of_freedom: number of additional parameters at each end; reasonable values depend on the kernel span +template +PiecewiseFittedCurve fit_curve( + const std::vector> &observations, + const std::vector &observation_points, + const std::vector &weights, + size_t segments_count, + size_t endpoints_level_of_freedom) { + + // check to make sure inputs are correct + assert(segments_count > 0); + assert(observations.size() > 0); + assert(observation_points.size() == observations.size()); + assert(observation_points.size() == weights.size()); + assert(segments_count <= observations.size()); + + //prepare sqrt of weights, which will then be applied to both matrix T and observed data: https://en.wikipedia.org/wiki/Weighted_least_squares + std::vector sqrt_weights(weights.size()); + for (size_t index = 0; index < weights.size(); ++index) { + assert(weights[index] > 0); + sqrt_weights[index] = sqrt(weights[index]); + } + + // prepare result and compute metadata + PiecewiseFittedCurve result { }; + + NumberType valid_length = observation_points.back() - observation_points.front(); + NumberType segment_size = valid_length / NumberType(segments_count); + result.start = observation_points.front(); + result.segment_size = segment_size; + result.endpoints_level_of_freedom = endpoints_level_of_freedom; + + // prepare observed data + // Eigen defaults to column major memory layout. + Eigen::MatrixXf data_points(Dimension, observations.size()); + for (size_t index = 0; index < observations.size(); ++index) { + data_points.col(index) = observations[index] * sqrt_weights[index]; + } + // parameters count is always increased by one to make the parametric space of the curve symmetric. + // without this fix, the end of the curve is less flexible than the beginning + size_t parameters_count = segments_count + 1 + 2 * endpoints_level_of_freedom; + //Create weight matrix T for each point and each segment; + Eigen::MatrixXf T(observation_points.size(), parameters_count); + T.setZero(); + //Fill the weight matrix + for (size_t i = 0; i < observation_points.size(); ++i) { + NumberType observation_point = observation_points[i]; + //find corresponding segment index; expects kernels to be centered + int middle_right_segment_index = floor((observation_point - result.start) / result.segment_size); + //find index of first segment that is affected by the point i; this can be deduced from kernel_span + int start_segment_idx = middle_right_segment_index - Kernel::kernel_span / 2 + 1; + for (int segment_index = start_segment_idx; segment_index < int(start_segment_idx + Kernel::kernel_span); + segment_index++) { + NumberType segment_start = result.start + segment_index * result.segment_size; + NumberType normalized_segment_distance = (segment_start - observation_point) / result.segment_size; + + int parameter_index = segment_index + endpoints_level_of_freedom; + parameter_index = std::clamp(parameter_index, 0, int(parameters_count) - 1); + T(i, parameter_index) += Kernel::kernel(normalized_segment_distance) * sqrt_weights[i]; + } + } + +#ifdef LSQR_DEBUG + std::cout << "weight matrix: " << std::endl; + for (int obs = 0; obs < observation_points.size(); ++obs) { + std::cout << std::endl; + for (int segment = 0; segment < parameters_count; ++segment) { + std::cout << T(obs, segment) << " "; + } + } + std::cout << std::endl; +#endif + + // Solve for linear least square fit + result.coefficients.resize(Dimension, parameters_count); + const auto QR = T.fullPivHouseholderQr(); + for (size_t dim = 0; dim < Dimension; ++dim) { + result.coefficients.row(dim) = QR.solve(data_points.row(dim).transpose()); + } + + return result; +} + +template +PiecewiseFittedCurve> +fit_cubic_bspline( + const std::vector> &observations, + std::vector observation_points, + std::vector weights, + size_t segments_count, + size_t endpoints_level_of_freedom = 0) { + return fit_curve>(observations, observation_points, weights, segments_count, + endpoints_level_of_freedom); +} + +template +PiecewiseFittedCurve> +fit_catmul_rom_spline( + const std::vector> &observations, + std::vector observation_points, + std::vector weights, + size_t segments_count, + size_t endpoints_level_of_freedom = 0) { + return fit_curve>(observations, observation_points, weights, segments_count, + endpoints_level_of_freedom); +} + +} +} + +#endif /* SRC_LIBSLIC3R_GEOMETRY_CURVES_HPP_ */ diff --git a/src/libslic3r/KDTreeIndirect.hpp b/src/libslic3r/KDTreeIndirect.hpp index 12e462569..37c10827b 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/QuadricEdgeCollapse.cpp b/src/libslic3r/QuadricEdgeCollapse.cpp index 48a72e6d8..be0a90ec1 100644 --- a/src/libslic3r/QuadricEdgeCollapse.cpp +++ b/src/libslic3r/QuadricEdgeCollapse.cpp @@ -595,8 +595,7 @@ bool QuadricEdgeCollapse::is_flipped(const Vec3f & new_vertex, const EdgeInfos & e_infos, const indexed_triangle_set &its) { - static const float thr_pos = 1.0f - std::numeric_limits::epsilon(); - static const float thr_neg = -thr_pos; + static const float triangle_beauty_threshold = 1.0f - std::numeric_limits::epsilon(); static const float dot_thr = 0.2f; // Value from simplify mesh cca 80 DEG // for each vertex triangles @@ -615,7 +614,16 @@ bool QuadricEdgeCollapse::is_flipped(const Vec3f & new_vertex, d2.normalize(); float dot = d1.dot(d2); - if (dot > thr_pos || dot < thr_neg) return true; + if (dot > triangle_beauty_threshold || dot < -triangle_beauty_threshold) { // OK, the new triangle is suspiciously ugly, but it can still be better than the original + const Vec3f &v_orig = its.vertices[t[(e_info.edge) % 3]]; + Vec3f d1_orig = vf - v_orig; + d1_orig.normalize(); + Vec3f d2_orig = vs - v_orig; + d2_orig.normalize(); + if (std::fabs(d1_orig.dot(d2_orig)) < std::fabs(dot)) { // original was not that ugly, so return flipped + return true; + } // else original triangle was worse than the new, so don't discard the new yet + } // IMPROVE: propagate new normal Vec3f n = d1.cross(d2); n.normalize(); diff --git a/src/libslic3r/SLA/bicubic.h b/src/libslic3r/SLA/bicubic.h deleted file mode 100644 index 870d00dbd..000000000 --- a/src/libslic3r/SLA/bicubic.h +++ /dev/null @@ -1,186 +0,0 @@ -#ifndef BICUBIC_HPP -#define BICUBIC_HPP - -#include -#include -#include - -#include - -namespace Slic3r { - -namespace BicubicInternal { - // Linear kernel, to be able to test cubic methods with hat kernels. - template - struct LinearKernel - { - typedef T FloatType; - - static T a00() { return T(0.); } - static T a01() { return T(0.); } - static T a02() { return T(0.); } - static T a03() { return T(0.); } - static T a10() { return T(1.); } - static T a11() { return T(-1.); } - static T a12() { return T(0.); } - static T a13() { return T(0.); } - static T a20() { return T(0.); } - static T a21() { return T(1.); } - static T a22() { return T(0.); } - static T a23() { return T(0.); } - static T a30() { return T(0.); } - static T a31() { return T(0.); } - static T a32() { return T(0.); } - static T a33() { return T(0.); } - }; - - // Interpolation kernel aka Catmul-Rom aka Keyes kernel. - template - struct CubicCatmulRomKernel - { - typedef T FloatType; - - static T a00() { return 0; } - static T a01() { return (T)-0.5; } - static T a02() { return (T) 1.; } - static T a03() { return (T)-0.5; } - static T a10() { return (T) 1.; } - static T a11() { return 0; } - static T a12() { return (T)-5./2.; } - static T a13() { return (T) 3./2.; } - static T a20() { return 0; } - static T a21() { return (T) 0.5; } - static T a22() { return (T) 2.; } - static T a23() { return (T)-3./2.; } - static T a30() { return 0; } - static T a31() { return 0; } - static T a32() { return (T)-0.5; } - static T a33() { return (T) 0.5; } - }; - - // B-spline kernel - template - struct CubicBSplineKernel - { - typedef T FloatType; - - static T a00() { return (T) 1./6.; } - static T a01() { return (T) -3./6.; } - static T a02() { return (T) 3./6.; } - static T a03() { return (T) -1./6.; } - static T a10() { return (T) 4./6.; } - static T a11() { return 0; } - static T a12() { return (T) -6./6.; } - static T a13() { return (T) 3./6.; } - static T a20() { return (T) 1./6.; } - static T a21() { return (T) 3./6.; } - static T a22() { return (T) 3./6.; } - static T a23() { return (T)- 3./6.; } - static T a30() { return 0; } - static T a31() { return 0; } - static T a32() { return 0; } - static T a33() { return (T) 1./6.; } - }; - - template - inline T clamp(T a, T lower, T upper) - { - return (a < lower) ? lower : - (a > upper) ? upper : a; - } -} - -template -struct CubicKernel -{ - typedef typename KERNEL KernelInternal; - typedef typename KERNEL::FloatType FloatType; - - static FloatType kernel(FloatType x) - { - x = fabs(x); - if (x >= (FloatType)2.) - return 0.0f; - if (x <= (FloatType)1.) { - FloatType x2 = x * x; - FloatType x3 = x2 * x; - return KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3; - } - assert(x > (FloatType)1. && x < (FloatType)2.); - x -= (FloatType)1.; - FloatType x2 = x * x; - FloatType x3 = x2 * x; - return KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3; - } - - static FloatType interpolate(FloatType f0, FloatType f1, FloatType f2, FloatType f3, FloatType x) - { - const FloatType x2 = x*x; - const FloatType x3 = x*x*x; - return f0*(KERNEL::a00() + KERNEL::a01() * x + KERNEL::a02() * x2 + KERNEL::a03() * x3) + - f1*(KERNEL::a10() + KERNEL::a11() * x + KERNEL::a12() * x2 + KERNEL::a13() * x3) + - f2*(KERNEL::a20() + KERNEL::a21() * x + KERNEL::a22() * x2 + KERNEL::a23() * x3) + - f3*(KERNEL::a30() + KERNEL::a31() * x + KERNEL::a32() * x2 + KERNEL::a33() * x3); - } -}; - -// Linear splines -typedef CubicKernel> LinearKernelf; -typedef CubicKernel> LinearKerneld; -// Catmul-Rom splines -typedef CubicKernel> CubicCatmulRomKernelf; -typedef CubicKernel> CubicCatmulRomKerneld; -typedef CubicKernel> CubicInterpolationKernelf; -typedef CubicKernel> CubicInterpolationKerneld; -// Cubic B-splines -typedef CubicKernel> CubicBSplineKernelf; -typedef CubicKernel> CubicBSplineKerneld; - -template -static float cubic_interpolate(const Eigen::ArrayBase &F, const typename KERNEL::FloatType pt, const typename KERNEL::FloatType dx) -{ - typedef typename KERNEL::FloatType T; - const int w = int(F.size()); - const int ix = (int)floor(pt); - const T s = pt - (T)ix; - - if (ix > 1 && ix + 2 < w) { - // Inside the fully interpolated region. - return KERNEL::interpolate(F[ix - 1], F[ix], F[ix + 1], F[ix + 2], s); - } - // Transition region. Extend with a constant function. - auto f = [&F, w](x) { return F[BicubicInternal::clamp(x, 0, w - 1)]; } - return KERNEL::interpolate(f(ix - 1), f(ix), f(ix + 1), f(ix + 2), s); -} - -template -static float bicubic_interpolate(const Eigen::MatrixBase &F, const Eigen::Matrix &pt, const typename KERNEL::FloatType dx) -{ - typedef typename KERNEL::FloatType T; - const int w = F.cols(); - const int h = F.rows(); - const int ix = (int)floor(pt[0]); - const int iy = (int)floor(pt[1]); - const T s = pt[0] - (T)ix; - const T t = pt[1] - (T)iy; - - if (ix > 1 && ix + 2 < w && iy > 1 && iy + 2 < h) { - // Inside the fully interpolated region. - return KERNEL::interpolate( - KERNEL::interpolate(F(ix-1,iy-1),F(ix ,iy-1),F(ix+1,iy-1),F(ix+2,iy-1),s), - KERNEL::interpolate(F(ix-1,iy ),F(ix ,iy ),F(ix+1,iy ),F(ix+2,iy ),s), - KERNEL::interpolate(F(ix-1,iy+1),F(ix ,iy+1),F(ix+1,iy+1),F(ix+2,iy+1),s), - KERNEL::interpolate(F(ix-1,iy+2),F(ix ,iy+2),F(ix+1,iy+2),F(ix+2,iy+2),s),t); - } - // Transition region. Extend with a constant function. - auto f = [&f, w, h](int x, int y) { return F(BicubicInternal::clamp(x,0,w-1),BicubicInternal::clamp(y,0,h-1)); } - return KERNEL::interpolate( - KERNEL::interpolate(f(ix-1,iy-1),f(ix ,iy-1),f(ix+1,iy-1),f(ix+2,iy-1),s), - KERNEL::interpolate(f(ix-1,iy ),f(ix ,iy ),f(ix+1,iy ),f(ix+2,iy ),s), - KERNEL::interpolate(f(ix-1,iy+1),f(ix ,iy+1),f(ix+1,iy+1),f(ix+2,iy+1),s), - KERNEL::interpolate(f(ix-1,iy+2),f(ix ,iy+2),f(ix+1,iy+2),f(ix+2,iy+2),s),t); -} - -} // namespace Slic3r - -#endif /* BICUBIC_HPP */ diff --git a/src/libslic3r/Subdivide.cpp b/src/libslic3r/Subdivide.cpp new file mode 100644 index 000000000..31988a851 --- /dev/null +++ b/src/libslic3r/Subdivide.cpp @@ -0,0 +1,218 @@ +#include "Subdivide.hpp" +#include "Point.hpp" + +namespace Slic3r{ + +indexed_triangle_set its_subdivide( + const indexed_triangle_set &its, float max_length) +{ + // same order as key order in Edge Divides + struct VerticesSequence + { + size_t start_index; + bool positive_order; + VerticesSequence(size_t start_index, bool positive_order = true) + : start_index(start_index), positive_order(positive_order){} + }; + // vertex index small, big vertex index from key.first to key.second + using EdgeDivides = std::map, VerticesSequence>; + struct Edges + { + Vec3f data[3]; + Vec3f lengths; + Edges(const Vec3crd &indices, const std::vector &vertices) + : lengths(-1.f,-1.f,-1.f) + { + const Vec3f &v0 = vertices[indices[0]]; + const Vec3f &v1 = vertices[indices[1]]; + const Vec3f &v2 = vertices[indices[2]]; + data[0] = v0 - v1; + data[1] = v1 - v2; + data[2] = v2 - v0; + } + float abs_sum(const Vec3f &v) + { + return abs(v[0]) + abs(v[1]) + abs(v[2]); + } + bool is_dividable(const float& max_length) { + Vec3f sum(abs_sum(data[0]), abs_sum(data[1]), abs_sum(data[2])); + Vec3i biggest_index = (sum[0] > sum[1]) ? + ((sum[0] > sum[2]) ? + ((sum[2] > sum[1]) ? + Vec3i(0, 2, 1) : + Vec3i(0, 1, 2)) : + Vec3i(2, 0, 1)) : + ((sum[1] > sum[2]) ? + ((sum[2] > sum[0]) ? + Vec3i(1, 2, 0) : + Vec3i(1, 0, 2)) : + Vec3i(2, 1, 0)); + for (int i = 0; i < 3; i++) { + int index = biggest_index[i]; + if (sum[index] <= max_length) return false; + lengths[index] = data[index].norm(); + if (lengths[index] <= max_length) continue; + + // calculate rest of lengths + for (int j = i + 1; j < 3; j++) { + index = biggest_index[j]; + lengths[index] = data[index].norm(); + } + return true; + } + return false; + } + }; + struct TriangleLengths + { + Vec3crd indices; + Vec3f l; // lengths + TriangleLengths(const Vec3crd &indices, const Vec3f &lengths) + : indices(indices), l(lengths) + {} + + int get_divide_index(float max_length) { + if (l[0] > l[1] && l[0] > l[2]) { + if (l[0] > max_length) return 0; + } else if (l[1] > l[2]) { + if (l[1] > max_length) return 1; + } else { + if (l[2] > max_length) return 2; + } + return -1; + } + + // divide triangle add new vertex to vertices + std::pair divide( + int divide_index, float max_length, + std::vector &vertices, + EdgeDivides &edge_divides) + { + // index to lengths and indices + size_t i0 = divide_index; + size_t i1 = (divide_index + 1) % 3; + size_t vi0 = indices[i0]; + size_t vi1 = indices[i1]; + std::pair key(vi0, vi1); + bool key_swap = false; + if (key.first > key.second) { + std::swap(key.first, key.second); + key_swap = true; + } + + float length = l[divide_index]; + size_t count_edge_vertices = static_cast(floor(length / max_length)); + float count_edge_segments = static_cast(count_edge_vertices + 1); + + auto it = edge_divides.find(key); + if (it == edge_divides.end()) { + // Create new vertices + VerticesSequence new_vs(vertices.size()); + Vec3f vf = vertices[key.first]; // copy + const Vec3f &vs = vertices[key.second]; + Vec3f dir = vs - vf; + for (size_t i = 1; i <= count_edge_vertices; ++i) { + float ratio = i / count_edge_segments; + vertices.push_back(vf + dir * ratio); + } + bool success; + std::tie(it,success) = edge_divides.insert({key, new_vs}); + assert(success); + } + const VerticesSequence &vs = it->second; + + int index_offset = count_edge_vertices/2; + size_t i2 = (divide_index + 2) % 3; + if (count_edge_vertices % 2 == 0 && key_swap == (l[i1] < l[i2])) { + --index_offset; + } + int sign = (vs.positive_order) ? 1 : -1; + size_t new_index = vs.start_index + sign*index_offset; + + size_t vi2 = indices[i2]; + const Vec3f &v2 = vertices[vi2]; + Vec3f new_edge = v2 - vertices[new_index]; + float new_len = new_edge.norm(); + + float ratio = (1 + index_offset) / count_edge_segments; + float len1 = l[i0] * ratio; + float len2 = l[i0] - len1; + if (key_swap) std::swap(len1, len2); + + Vec3crd indices1(vi0, new_index, vi2); + Vec3f lengths1(len1, new_len, l[i2]); + + Vec3crd indices2(new_index, vi1, vi2); + Vec3f lengths2(len2, l[i1], new_len); + + // append key for divided edge when neccesary + if (index_offset > 0) { + std::pair new_key(key.first, new_index); + bool new_key_swap = false; + if (new_key.first > new_key.second) { + std::swap(new_key.first, new_key.second); + new_key_swap = true; + } + if (edge_divides.find(new_key) == edge_divides.end()) { + // insert new + edge_divides.insert({new_key, (new_key_swap) ? + VerticesSequence(new_index - sign, !vs.positive_order) + : vs}); + } + } + + if (index_offset < int(count_edge_vertices)-1) { + std::pair new_key(new_index, key.second); + bool new_key_swap = false; + if (new_key.first > new_key.second) { + std::swap(new_key.first, new_key.second); + new_key_swap = true; + } + // bad order + if (edge_divides.find(new_key) == edge_divides.end()) { + edge_divides.insert({new_key, (new_key_swap) ? + VerticesSequence(vs.start_index + sign*(count_edge_vertices-1), !vs.positive_order) + : VerticesSequence(new_index + sign, vs.positive_order)}); + } + } + + return {TriangleLengths(indices1, lengths1), + TriangleLengths(indices2, lengths2)}; + } + }; + indexed_triangle_set result; + result.indices.reserve(its.indices.size()); + const std::vector &vertices = its.vertices; + result.vertices = vertices; // copy + std::queue tls; + + EdgeDivides edge_divides; + for (const Vec3crd &indices : its.indices) { + Edges edges(indices, vertices); + // speed up only sum not sqrt is apply + if (!edges.is_dividable(max_length)) { + // small triangle + result.indices.push_back(indices); + continue; + } + TriangleLengths tl(indices, edges.lengths); + do { + int divide_index = tl.get_divide_index(max_length); + if (divide_index < 0) { + // no dividing + result.indices.push_back(tl.indices); + if (tls.empty()) break; + tl = tls.front(); // copy + tls.pop(); + } else { + auto [tl1, tl2] = tl.divide(divide_index, max_length, + result.vertices, edge_divides); + tl = tl1; + tls.push(tl2); + } + } while (true); + } + return result; +} + +} diff --git a/src/libslic3r/Subdivide.hpp b/src/libslic3r/Subdivide.hpp new file mode 100644 index 000000000..f97e4b3c5 --- /dev/null +++ b/src/libslic3r/Subdivide.hpp @@ -0,0 +1,12 @@ +#ifndef libslic3r_Subdivide_hpp_ +#define libslic3r_Subdivide_hpp_ + +#include "TriangleMesh.hpp" + +namespace Slic3r { + +indexed_triangle_set its_subdivide(const indexed_triangle_set &its, float max_length); + +} + +#endif //libslic3r_Subdivide_hpp_ diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 5e0fb67b3..06cad840e 100644 --- a/src/libslic3r/libslic3r.h +++ b/src/libslic3r/libslic3r.h @@ -23,6 +23,13 @@ #include #include +#ifdef _WIN32 +// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation. +// https://github.com/microsoft/STL/issues/147#issuecomment-1090148740 +// Thus it is recommended to use boost::container::deque instead. +#include +#endif // _WIN32 + #include "Technologies.hpp" #include "Semver.hpp" @@ -73,6 +80,16 @@ namespace Slic3r { extern Semver SEMVER; +// On MSVC, std::deque degenerates to a list of pointers, which defeats its purpose of reducing allocator load and memory fragmentation. +template> +using deque = +#ifdef _WIN32 + // Use boost implementation, which allocates blocks of 512 bytes instead of blocks of 8 bytes. + boost::container::deque; +#else // _WIN32 + std::deque; +#endif // _WIN32 + template inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); } @@ -314,6 +331,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 From 6da220062ca4e07ef670ad6495d5aad8864d7bf8 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Mon, 23 May 2022 10:30:45 +0200 Subject: [PATCH 02/11] Mac OS pre 10.13 does not fully support std::optional (method .value() is not allowed) This commit replaces usage of .value() calls with .operator*() --- src/libslic3r/GCode/SeamPlacer.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index 4f30fdc19..d84fc00f7 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -1179,17 +1179,18 @@ std::vector> SeamPlacer::find_seam_string(const PrintO float(po->get_layer(next_layer)->slice_z), comparator); if (maybe_next_seam.has_value()) { - std::pair next_seam_coords = maybe_next_seam.value(); + // For old macOS (pre 10.14), std::optional does not have .value() method, so the code is using operator*() instead. + std::pair next_seam_coords = maybe_next_seam.operator*(); const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second]; bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second; out_moved_seams_count += is_moved; if (is_moved && (!out_best_moved_seam.has_value() || comparator.is_first_better(next_seam, - layers[out_best_moved_seam.value().first].points[out_best_moved_seam.value().second]))) { + layers[out_best_moved_seam.operator*().first].points[out_best_moved_seam.operator*().second]))) { out_best_moved_seam = { next_seam_coords }; } - seam_string.push_back(maybe_next_seam.value()); + seam_string.push_back(maybe_next_seam.operator*()); prev_point_index = seam_string.back(); //String added, prev_point_index updated } else { @@ -1208,17 +1209,17 @@ std::vector> SeamPlacer::find_seam_string(const PrintO float(po->get_layer(next_layer)->slice_z), comparator); if (maybe_next_seam.has_value()) { - std::pair next_seam_coords = maybe_next_seam.value(); + std::pair next_seam_coords = maybe_next_seam.operator*(); const auto &next_seam = layers[next_seam_coords.first].points[next_seam_coords.second]; bool is_moved = next_seam.perimeter.seam_index != next_seam_coords.second; out_moved_seams_count += is_moved; if (is_moved && (!out_best_moved_seam.has_value() || comparator.is_first_better(next_seam, - layers[out_best_moved_seam.value().first].points[out_best_moved_seam.value().second]))) { + layers[out_best_moved_seam.operator*().first].points[out_best_moved_seam.operator*().second]))) { out_best_moved_seam = { next_seam_coords }; } - seam_string.push_back(maybe_next_seam.value()); + seam_string.push_back(maybe_next_seam.operator*()); prev_point_index = seam_string.back(); //String added, prev_point_index updated } else { @@ -1303,7 +1304,7 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl:: moved_seams_count); if (best_moved_seam.has_value()) { size_t alternative_moved_seams_count; - alternative_seam_string = this->find_seam_string(po, best_moved_seam.value(), comparator, + alternative_seam_string = this->find_seam_string(po, best_moved_seam.operator*(), comparator, best_moved_seam, alternative_moved_seams_count); if (alternative_seam_string.size() >= SeamPlacer::seam_align_minimum_string_seams && alternative_moved_seams_count < moved_seams_count) { From a4201321e8c66fdbd27d6c8ede42794fcac03610 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Wed, 25 May 2022 16:56:45 +0200 Subject: [PATCH 03/11] Hopefully improved the seam placer performance a lot --- src/libslic3r/CMakeLists.txt | 2 + src/libslic3r/GCode/SeamPlacer.cpp | 286 +++++++++++++------------- src/libslic3r/GCode/SeamPlacer.hpp | 20 +- src/libslic3r/TriangleSetSampling.cpp | 72 +++++++ src/libslic3r/TriangleSetSampling.hpp | 20 ++ 5 files changed, 243 insertions(+), 157 deletions(-) create mode 100644 src/libslic3r/TriangleSetSampling.cpp create mode 100644 src/libslic3r/TriangleSetSampling.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 21fac7456..4f7044379 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -245,6 +245,8 @@ add_library(libslic3r STATIC Thread.hpp TriangleSelector.cpp TriangleSelector.hpp + TriangleSetSampling.cpp + TriangleSetSampling.hpp MTUtils.hpp Zipper.hpp Zipper.cpp diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index d84fc00f7..a9f8dd803 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -9,15 +9,15 @@ #include #include "libslic3r/AABBTreeLines.hpp" +#include "libslic3r/KDTreeIndirect.hpp" #include "libslic3r/ExtrusionEntity.hpp" #include "libslic3r/Print.hpp" #include "libslic3r/BoundingBox.hpp" #include "libslic3r/ClipperUtils.hpp" #include "libslic3r/Layer.hpp" -#include "libslic3r/QuadricEdgeCollapse.hpp" -#include "libslic3r/Subdivide.hpp" #include "libslic3r/Geometry/Curves.hpp" +#include "libslic3r/TriangleSetSampling.hpp" #include "libslic3r/Utils.hpp" @@ -130,19 +130,22 @@ Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) { return Vec3f(cos(term1) * term3, sin(term1) * term3, term2); } -std::vector raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree, - const indexed_triangle_set &triangles, size_t negative_volumes_start_index) { +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 for " << triangles.indices.size() << " triangles: start"; + << "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_triangle; + float step_size = 1.0f / SeamPlacer::sqr_rays_per_sample_point; std::vector precomputed_sample_directions( - SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); - for (size_t x_idx = 0; x_idx < SeamPlacer::sqr_rays_per_triangle; ++x_idx) { + 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_triangle; ++y_idx) { - size_t dir_index = x_idx * SeamPlacer::sqr_rays_per_triangle + y_idx; + 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 }); } @@ -150,24 +153,19 @@ std::vector raycast_visibility(const AABBTreeIndirect::Tree< bool model_contains_negative_parts = negative_volumes_start_index < triangles.indices.size(); - std::vector result(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](tbb::blocked_range r) { + &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 face_index = r.begin(); face_index < r.end(); ++face_index) { - FaceVisibilityInfo &dest = result[face_index]; - dest.visibility = 1.0f; - constexpr float decrease = 1.0f - / (SeamPlacer::sqr_rays_per_triangle * SeamPlacer::sqr_rays_per_triangle); + 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); - Vec3i face = triangles.indices[face_index]; - Vec3f A = triangles.vertices[face.x()]; - Vec3f B = triangles.vertices[face.y()]; - Vec3f C = triangles.vertices[face.z()]; - Vec3f center = (A + B + C) / 3.0f; - Vec3f normal = ((B - A).cross(C - B)).normalized(); + 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); @@ -183,11 +181,14 @@ std::vector raycast_visibility(const AABBTreeIndirect::Tree< 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) { - dest.visibility -= decrease; + 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; + Vec3d ray_origin_d = (center + normal * 0.1).cast(); // start above surface. - if (face_index >= negative_volumes_start_index) { // if casting from negative volume face, invert direction, change start pos + 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.03).cast(); } @@ -209,7 +210,7 @@ std::vector raycast_visibility(const AABBTreeIndirect::Tree< } } if (counter == 0) { - dest.visibility -= decrease; + result[s_idx] -= decrease_step; } } } @@ -218,7 +219,8 @@ std::vector raycast_visibility(const AABBTreeIndirect::Tree< }); BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end"; + << "SeamPlacer: raycast visibility of " << samples.positions.size() << " samples over " << triangles.indices.size() + << " triangles: end"; return result; } @@ -273,12 +275,28 @@ std::vector calculate_polygon_angles_at_vertices(const Polygon &polygon, return result; } +struct CoordinateFunctor { + const std::vector *coordinates; + CoordinateFunctor(const std::vector *coords) : + coordinates(coords) { + } + CoordinateFunctor() : + coordinates(nullptr) { + } + + const float& operator()(size_t idx, size_t dim) const { + return coordinates->operator [](idx)[dim]; + } +}; + // structure to store global information about the model - occlusion hits, enforcers, blockers struct GlobalModelInfo { - indexed_triangle_set model; - std::vector triangle_neighbours; - AABBTreeIndirect::Tree<3, float> model_tree; - std::vector visiblity_info; + 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; + indexed_triangle_set enforcers; indexed_triangle_set blockers; AABBTreeIndirect::Tree<3, float> enforcers_tree; @@ -303,48 +321,72 @@ struct GlobalModelInfo { } float calculate_point_visibility(const Vec3f &position) const { - size_t hit_idx; - Vec3f hit_point; - if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree, - position, hit_idx, hit_point) >= 0) { - float visibility = visiblity_info[hit_idx].visibility; - Vec3i neighbours = this->triangle_neighbours[hit_idx]; - size_t n_count = 0; - for (int neighbour : neighbours) { - if (neighbour >= 0) { - visibility += visiblity_info[neighbour].visibility; - n_count++; - } - } - return visibility / (1 + n_count); - } else { - return 0.0f; + std::vector points = find_nearby_points(mesh_samples_tree, position, mesh_samples_radius); + if (points.empty()) { + size_t idx = find_closest_point(mesh_samples_tree, position); + return mesh_samples_visibility[idx]; } + + float total_weight = 0; + float total_visibility = 0; + for (size_t i = 0; i < points.size(); ++i) { + size_t sample_idx = points[i]; + float weight = 1.0f; // SeamPlacer::visibility_samples_radius * SeamPlacer::visibility_samples_radius - + //(position - mesh_samples.positions[sample_idx]).squaredNorm(); + total_visibility += weight * mesh_samples_visibility[sample_idx]; + total_weight += weight; + } + + return total_visibility / total_weight; + } #ifdef DEBUG_FILES - void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const { + void debug_export(const indexed_triangle_set &obj_mesh) const { + indexed_triangle_set divided_mesh = obj_mesh; Slic3r::CNumericLocalesSetter locales_setter; - FILE *fp = boost::nowide::fopen(file_name, "w"); - if (fp == nullptr) { - BOOST_LOG_TRIVIAL(error) - << "stl_write_obj: Couldn't open " << file_name << " for writing"; - return; + { + 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); } - 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)); + { + 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); } - for (size_t i = 0; i < divided_mesh.indices.size(); ++i) - fprintf(fp, "f %d %d %d\n", divided_mesh.indices[i][0] + 1, divided_mesh.indices[i][1] + 1, - divided_mesh.indices[i][2] + 1); - fclose(fp); } #endif @@ -563,86 +605,40 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { } } } - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: gather occlusion meshes: end"; - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: simplify occlusion meshes: start"; - - //simplify raycasting mesh - { - its_quadric_edge_collapse(triangle_set, SeamPlacer::raycasting_decimation_target_triangle_count, nullptr, - nullptr, - nullptr); - float triangle_set_area = tbb::parallel_reduce(tbb::blocked_range(0, triangle_set.indices.size()), 0, - [&triangle_set]( - tbb::blocked_range r, float sum) { - 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()]; - sum += 0.5f * (b - a).cross(c - a).norm(); - } - return sum; - }, std::plus()); - - float target_triangle_area = triangle_set_area / SeamPlacer::raycasting_subdivision_target_triangle_count; - float target_triangle_length = 2 * 1.316 * sqrtf(target_triangle_area); //assuming 30-30-120 triangle - float subdivision_length = std::max(SeamPlacer::raycasting_subdivision_target_length, target_triangle_length); - triangle_set = its_subdivide(triangle_set, subdivision_length); - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: triangle set after subdivision: " << triangle_set.indices.size(); - } - - //simplify negative volumes - { - its_quadric_edge_collapse(negative_volumes_set, SeamPlacer::raycasting_decimation_target_triangle_count, - nullptr, - nullptr, - nullptr); - float negative_volumes_set_area = tbb::parallel_reduce( - tbb::blocked_range(0, negative_volumes_set.indices.size()), 0, - [&negative_volumes_set]( - tbb::blocked_range r, float sum) { - for (size_t t_idx = r.begin(); t_idx < r.end(); ++t_idx) { - const Vec3f &a = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].x()]; - const Vec3f &b = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].y()]; - const Vec3f &c = negative_volumes_set.vertices[negative_volumes_set.indices[t_idx].z()]; - sum += 0.5f * (b - a).cross(c - a).norm(); - } - return sum; - }, std::plus()); - - float target_triangle_area = negative_volumes_set_area - / SeamPlacer::raycasting_subdivision_target_triangle_count; - float target_triangle_length = 2 * 1.316 * sqrtf(target_triangle_area); //assuming 30-30-120 triangle - float subdivision_length = std::max(SeamPlacer::raycasting_subdivision_target_length, target_triangle_length); - negative_volumes_set = its_subdivide(negative_volumes_set, subdivision_length); - } size_t negative_volumes_start_index = triangle_set.indices.size(); its_merge(triangle_set, negative_volumes_set); its_transform(triangle_set, obj_transform); BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: simplify occlusion meshes: end"; + << "SeamPlacer: gather occlusion meshes: end"; + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: Compute visiblity 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()); + result.mesh_samples_radius = sqrt( + 4.0f * (result.mesh_samples.total_area / SeamPlacer::raycasting_visibility_samples_count) / PI); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: Compute visiblity sample points: end; mesh_sample_radius: " << 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); - std::vector neighbours = its_face_neighbors_par(triangle_set); BOOST_LOG_TRIVIAL(debug) << "SeamPlacer:build AABB tree: end"; - result.model = triangle_set; - result.triangle_neighbours = neighbours; - result.model_tree = raycasting_tree; - result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set, negative_volumes_start_index); + result.mesh_samples_visibility = raycast_visibility(raycasting_tree, triangle_set, result.mesh_samples, + negative_volumes_start_index); #ifdef DEBUG_FILES - auto filename = debug_out_path("visiblity.obj"); - result.debug_export(triangle_set, filename.c_str()); + result.debug_export(triangle_set); #endif } @@ -1408,26 +1404,28 @@ void SeamPlacer::init(const Print &print) { 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); + { + GlobalModelInfo global_model_info { }; + gather_enforcers_blockers(global_model_info, po); - if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { - compute_global_occlusion(global_model_info, po); - } + if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { + compute_global_occlusion(global_model_info, po); + } - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: gather_seam_candidates: start"; - gather_seam_candidates(po, global_model_info, configured_seam_preference); - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: gather_seam_candidates: end"; - - if (configured_seam_preference == spAligned || configured_seam_preference == spNearest) { BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: calculate_candidates_visibility : start"; - calculate_candidates_visibility(po, global_model_info); + << "SeamPlacer: gather_seam_candidates: start"; + gather_seam_candidates(po, global_model_info, configured_seam_preference); BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: calculate_candidates_visibility : end"; - } + << "SeamPlacer: gather_seam_candidates: end"; + + 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) BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: calculate_overhangs and layer embdedding : start"; @@ -1495,7 +1493,7 @@ void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool extern Vec3f seam_position; size_t seam_index; - if (const Perimeter &perimeter = layer_perimeters.points[closest_perimeter_point_index].perimeter ; + 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; diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 7dbe7bcb3..c432668a3 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -87,10 +87,6 @@ struct SeamCandidate { bool central_enforcer; //marks this candidate as central point of enforced segment on the perimeter - important for alignment }; -struct FaceVisibilityInfo { - float visibility; -}; - struct SeamCandidateCoordinateFunctor { SeamCandidateCoordinateFunctor(const std::vector &seam_candidates) : seam_candidates(seam_candidates) { @@ -125,15 +121,13 @@ struct PrintObjectSeamData class SeamPlacer { public: - static constexpr size_t raycasting_decimation_target_triangle_count = 10000; - // for subdivision, both following criteria are considered, and the one with less resulting triangles is used - static constexpr size_t raycasting_subdivision_target_triangle_count = 20000; - static constexpr float raycasting_subdivision_target_length = 2.0f; - //square of number of rays per triangle - static constexpr size_t sqr_rays_per_triangle = 7; + // 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 = 40000; + //square of number of rays per sample point + static constexpr size_t sqr_rays_per_sample_point = 8; // arm length used during angles computation - static constexpr float polygon_local_angles_arm_distance = 0.1f; + static constexpr float polygon_local_angles_arm_distance = 0.3f; // max tolerable distance from the previous layer is overhang_distance_tolerance_factor * flow_width @@ -141,7 +135,7 @@ public: // determines angle importance compared to visibility ( neutral value is 1.0f. ) - static constexpr float angle_importance = 0.7f; + static constexpr float angle_importance = 0.6f; // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer static constexpr float enforcer_blocker_distance_tolerance = 0.35f; @@ -150,7 +144,7 @@ public: // 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.25f; + static constexpr float seam_align_score_tolerance = 0.27f; // seam_align_tolerable_dist - if next layer closes point is too far away, break string static constexpr float seam_align_tolerable_dist = 1.0f; // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. diff --git a/src/libslic3r/TriangleSetSampling.cpp b/src/libslic3r/TriangleSetSampling.cpp new file mode 100644 index 000000000..36022d8e3 --- /dev/null +++ b/src/libslic3r/TriangleSetSampling.cpp @@ -0,0 +1,72 @@ +#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()]; + float area = 0.5f * (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::random_device rnd_device; + std::mt19937 mersenne_engine { rnd_device() }; + // random numbers on interval [0, 1) + std::uniform_real_distribution fdistribution; + + auto get_random = [&fdistribution, &mersenne_engine]() { + return Vec3f { 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) { + float t_sample = random_samples[s_idx].x() * area_sum; + size_t t_idx = area_sum_to_triangle_idx.upper_bound(t_sample)->second; + + float sq_u = std::sqrt(random_samples[s_idx].y()); + float 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 000000000..28a661d76 --- /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_ */ From c23d1488c9b41be7a10a4b592d44ffd7301f4d66 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Fri, 27 May 2022 13:15:29 +0200 Subject: [PATCH 04/11] Performance improvements --- src/libslic3r/GCode/SeamPlacer.cpp | 15 +++++++++------ src/libslic3r/GCode/SeamPlacer.hpp | 6 +++--- src/libslic3r/TriangleSetSampling.cpp | 21 ++++++++++----------- 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index a9f8dd803..6def921de 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -323,8 +323,7 @@ struct GlobalModelInfo { float calculate_point_visibility(const Vec3f &position) const { std::vector points = find_nearby_points(mesh_samples_tree, position, mesh_samples_radius); if (points.empty()) { - size_t idx = find_closest_point(mesh_samples_tree, position); - return mesh_samples_visibility[idx]; + return 1.0f; } float total_weight = 0; @@ -622,18 +621,22 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { result.mesh_samples_tree = KDTreeIndirect<3, float, CoordinateFunctor>(result.mesh_samples_coordinate_functor, result.mesh_samples.positions.size()); result.mesh_samples_radius = sqrt( - 4.0f * (result.mesh_samples.total_area / SeamPlacer::raycasting_visibility_samples_count) / PI); + 10.0f * (result.mesh_samples.total_area / SeamPlacer::raycasting_visibility_samples_count) / PI); BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: Compute visiblity sample points: end; mesh_sample_radius: " << result.mesh_samples_radius; + << "SeamPlacer: Compute visiblity sample points: end"; + BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer:build AABB tree: start"; + << "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); BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer:build AABB tree: end"; + << "SeamPlacer: build AABB tree: end"; result.mesh_samples_visibility = raycast_visibility(raycasting_tree, triangle_set, result.mesh_samples, negative_volumes_start_index); diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index c432668a3..0a38e9a66 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -122,7 +122,7 @@ struct PrintObjectSeamData 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 = 40000; + static constexpr size_t raycasting_visibility_samples_count = 25000; //square of number of rays per sample point static constexpr size_t sqr_rays_per_sample_point = 8; @@ -135,7 +135,7 @@ public: // determines angle importance compared to visibility ( neutral value is 1.0f. ) - static constexpr float angle_importance = 0.6f; + static constexpr float angle_importance = 0.5f; // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer static constexpr float enforcer_blocker_distance_tolerance = 0.35f; @@ -144,7 +144,7 @@ public: // 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.27f; + static constexpr float seam_align_score_tolerance = 0.3f; // seam_align_tolerable_dist - if next layer closes point is too far away, break string static constexpr float seam_align_tolerable_dist = 1.0f; // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. diff --git a/src/libslic3r/TriangleSetSampling.cpp b/src/libslic3r/TriangleSetSampling.cpp index 36022d8e3..30f9a9b40 100644 --- a/src/libslic3r/TriangleSetSampling.cpp +++ b/src/libslic3r/TriangleSetSampling.cpp @@ -7,7 +7,7 @@ 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()); + std::vector triangles_area(triangle_set.indices.size()); tbb::parallel_for(tbb::blocked_range(0, triangle_set.indices.size()), [&triangle_set, &triangles_area]( @@ -16,28 +16,27 @@ TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const index 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()]; - float area = 0.5f * (b - a).cross(c - a).norm(); + double area = double(0.5 * (b - a).cross(c - a).norm()); triangles_area[t_idx] = area; } }); - std::map area_sum_to_triangle_idx; + 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::random_device rnd_device; - std::mt19937 mersenne_engine { rnd_device() }; + std::mt19937_64 mersenne_engine { }; // random numbers on interval [0, 1) - std::uniform_real_distribution fdistribution; + std::uniform_real_distribution fdistribution; auto get_random = [&fdistribution, &mersenne_engine]() { - return Vec3f { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) }; + return Vec3d { fdistribution(mersenne_engine), fdistribution(mersenne_engine), fdistribution(mersenne_engine) }; }; - std::vector random_samples(samples_count); + std::vector random_samples(samples_count); std::generate(random_samples.begin(), random_samples.end(), get_random); TriangleSetSamples result; @@ -50,11 +49,11 @@ TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const index [&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) { - float t_sample = random_samples[s_idx].x() * area_sum; + double t_sample = random_samples[s_idx].x() * area_sum; size_t t_idx = area_sum_to_triangle_idx.upper_bound(t_sample)->second; - float sq_u = std::sqrt(random_samples[s_idx].y()); - float v = random_samples[s_idx].z(); + 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()]; From b5b39195f4caeb468572906b859d7e3490454a7b Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Fri, 27 May 2022 14:16:55 +0200 Subject: [PATCH 05/11] Added throw_if_canceled callback to all slower sections --- src/libslic3r/GCode.cpp | 3 ++- src/libslic3r/GCode/SeamPlacer.cpp | 37 +++++++++++++++++++----------- src/libslic3r/GCode/SeamPlacer.hpp | 8 +++---- 3 files changed, 29 insertions(+), 19 deletions(-) diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index d91078164..cb14668fb 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -1347,7 +1347,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); if (! (has_wipe_tower && print.config().single_extruder_multi_material_priming)) { // Set initial extruder only after custom start G-code. diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index 6def921de..b28ea92ba 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -17,6 +17,7 @@ #include "libslic3r/Layer.hpp" #include "libslic3r/Geometry/Curves.hpp" +#include "libslic3r/QuadricEdgeCollapse.hpp" #include "libslic3r/TriangleSetSampling.hpp" #include "libslic3r/Utils.hpp" @@ -584,7 +585,7 @@ std::pair find_previous_and_next_perimeter_point(const std::vect } // Computes all global model info - transforms object, performs raycasting -void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { +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(); @@ -604,6 +605,7 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { } } } + throw_if_canceled(); size_t negative_volumes_start_index = triangle_set.indices.size(); its_merge(triangle_set, negative_volumes_set); @@ -613,7 +615,15 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { << "SeamPlacer: gather occlusion meshes: end"; BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: Compute visiblity sample points: start"; + << "SeamPlacer: decimate: start"; + float target_error = 1.0f; + its_quadric_edge_collapse(triangle_set, 0, &target_error, throw_if_canceled); + 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); @@ -625,7 +635,7 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { 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; @@ -635,11 +645,12 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { 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 @@ -1398,29 +1409,29 @@ void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl:: } -void SeamPlacer::init(const Print &print) { +void SeamPlacer::init(const Print &print, std::function throw_if_canceled_func) { 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); + 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"; @@ -1429,13 +1440,13 @@ void SeamPlacer::init(const Print &print) { << "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"; @@ -1456,7 +1467,7 @@ void SeamPlacer::init(const Print &print) { BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: pick_seam_point : end"; } - + throw_if_canceled_func(); if (configured_seam_preference == spAligned) { BOOST_LOG_TRIVIAL(debug) << "SeamPlacer: align_seam_points : start"; diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 0a38e9a66..de7e2f67c 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -129,13 +129,11 @@ public: // arm length used during angles computation static constexpr float polygon_local_angles_arm_distance = 0.3f; - // max tolerable distance from the previous layer is overhang_distance_tolerance_factor * flow_width static constexpr float overhang_distance_tolerance_factor = 0.5f; - // determines angle importance compared to visibility ( neutral value is 1.0f. ) - static constexpr float angle_importance = 0.5f; + static constexpr float angle_importance = 0.6f; // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer static constexpr float enforcer_blocker_distance_tolerance = 0.35f; @@ -144,7 +142,7 @@ public: // 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; + static constexpr float seam_align_score_tolerance = 0.27f; // seam_align_tolerable_dist - if next layer closes point is too far away, break string static constexpr float seam_align_tolerable_dist = 1.0f; // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. @@ -158,7 +156,7 @@ public: //The following data structures hold all perimeter points for all PrintObject. std::unordered_map m_seam_per_object; - void init(const Print &print); + 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; From 1e7b4a672071efeb069436387cb3825a85efa877 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Thu, 2 Jun 2022 13:53:40 +0200 Subject: [PATCH 06/11] Implementation of ShortEdgeCollapse Replaced QEC by edge collapse in occlusion computation --- src/libslic3r/CMakeLists.txt | 2 + src/libslic3r/GCode/SeamPlacer.cpp | 22 ++-- src/libslic3r/GCode/SeamPlacer.hpp | 4 +- src/libslic3r/ShortEdgeCollapse.cpp | 166 ++++++++++++++++++++++++++++ src/libslic3r/ShortEdgeCollapse.hpp | 16 +++ 5 files changed, 197 insertions(+), 13 deletions(-) create mode 100644 src/libslic3r/ShortEdgeCollapse.cpp create mode 100644 src/libslic3r/ShortEdgeCollapse.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 4f7044379..e41d1e3cc 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -205,6 +205,8 @@ add_library(libslic3r STATIC QuadricEdgeCollapse.cpp QuadricEdgeCollapse.hpp Semver.cpp + ShortEdgeCollapse.cpp + ShortEdgeCollapse.hpp ShortestPath.cpp ShortestPath.hpp SLAPrint.cpp diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index b28ea92ba..2a7448d8f 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -17,7 +17,7 @@ #include "libslic3r/Layer.hpp" #include "libslic3r/Geometry/Curves.hpp" -#include "libslic3r/QuadricEdgeCollapse.hpp" +#include "libslic3r/ShortEdgeCollapse.hpp" #include "libslic3r/TriangleSetSampling.hpp" #include "libslic3r/Utils.hpp" @@ -585,7 +585,8 @@ std::pair find_previous_and_next_perimeter_point(const std::vect } // Computes all global model info - transforms object, performs raycasting -void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po, std::function throw_if_canceled) { +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(); @@ -607,20 +608,19 @@ void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po, st } throw_if_canceled(); - 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: gather occlusion meshes: end"; BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: decimate: start"; - float target_error = 1.0f; - its_quadric_edge_collapse(triangle_set, 0, &target_error, throw_if_canceled); - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: decimate: end"; + << "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"; diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index de7e2f67c..08974bf5f 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -122,9 +122,9 @@ struct PrintObjectSeamData 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 = 25000; + 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 = 8; + 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; diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp new file mode 100644 index 000000000..049f58b8d --- /dev/null +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -0,0 +1,166 @@ +#include "ShortEdgeCollapse.hpp" + +#include +#include +#include +#include + +namespace Slic3r { + +void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) { + std::unordered_map vertices_index_mapping; // this map has two usages: + // in the first step, it contains mapping from original vertex index to final vertex index (which is different from original only for removed vertices) + + std::unordered_set vertices_to_remove; + std::unordered_set faces_to_remove; + std::vector triangles_neighbors; + std::vector face_normals; + std::vector vertex_normals; + std::vector face_indices; + std::random_device rd; + std::mt19937 generator(rd()); + + float decimation_ratio = 1.0f; + float edge_size = 0.2f; + size_t triangle_count = mesh.indices.size(); + + while (triangle_count > target_triangle_count) { + if (decimation_ratio < 0.4) { + edge_size *= 1.5f; + } + if (decimation_ratio < 0.05) { + edge_size *= 1.5f; + } + + float max_edge_len_squared = edge_size * edge_size; + triangles_neighbors = its_face_neighbors_par(mesh); + face_normals = its_face_normals(mesh); + + vertex_normals.resize(mesh.vertices.size()); + face_indices.resize(mesh.indices.size()); + 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]; + vertex_normals[t[0]] = n; + vertex_normals[t[1]] = n; + vertex_normals[t[2]] = n; + + face_indices[face_idx] = face_idx; + } + + std::vector min_vertex_score(mesh.vertices.size(), 1); + 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_score[t[0]] = std::min(min_vertex_score[t[0]], n.dot(vertex_normals[t[0]])); + min_vertex_score[t[1]] = std::min(min_vertex_score[t[1]], n.dot(vertex_normals[t[1]])); + min_vertex_score[t[2]] = std::min(min_vertex_score[t[2]], n.dot(vertex_normals[t[2]])); + } + + for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { + vertices_index_mapping[vertex_index] = vertex_index; + } + + std::shuffle(face_indices.begin(), face_indices.end(), generator); + + for (const size_t &face_idx : face_indices) { + if (faces_to_remove.find(face_idx) != faces_to_remove.end()) { + continue; + } + for (size_t edge_idx = 0; edge_idx < 3; ++edge_idx) { + size_t vertex_index_keep = mesh.indices[face_idx][edge_idx]; + size_t vertex_index_remove = mesh.indices[face_idx][(edge_idx + 1) % 3]; + + if ((mesh.vertices[vertex_index_keep] - mesh.vertices[vertex_index_remove]).squaredNorm() + > max_edge_len_squared) { + continue; + } + + if (min_vertex_score[vertex_index_remove] < min_vertex_score[vertex_index_keep]) { + size_t tmp = vertex_index_keep; + vertex_index_keep = vertex_index_remove; + vertex_index_remove = tmp; + } + + if (min_vertex_score[vertex_index_remove] < -1.0f) { + continue; + } + + if (vertices_to_remove.find(vertex_index_keep) != vertices_to_remove.end() + || vertices_to_remove.find(vertex_index_remove) != vertices_to_remove.end()) { + break; + } + + int neighbor_face_idx = triangles_neighbors[face_idx][edge_idx]; + if (neighbor_face_idx > 0 && faces_to_remove.find(neighbor_face_idx) != faces_to_remove.end()) { + continue; + } + + faces_to_remove.insert(face_idx); + faces_to_remove.insert(neighbor_face_idx); + vertices_to_remove.insert(vertex_index_remove); + vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep]; + min_vertex_score[vertex_index_keep] = -2.0; + break; + } + + } + + //flatten the mapping + for (auto &pair : vertices_index_mapping) { + while (vertices_index_mapping[pair.second] != pair.second) { + pair.second = vertices_index_mapping[pair.second]; + } + } + + std::vector new_vertices; + new_vertices.reserve(mesh.vertices.size() - vertices_to_remove.size()); + for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { + if (vertices_to_remove.find(vertex_index) == vertices_to_remove.end()) { //not removed + new_vertices.push_back(mesh.vertices[vertex_index]); + assert(vertices_index_mapping[vertex_index] == vertex_index); + vertices_index_mapping[vertex_index] = new_vertices.size() - 1; + } + } + + std::vector new_indices; + new_indices.reserve(mesh.indices.size() - faces_to_remove.size()); + for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) { + if (faces_to_remove.find(face_idx) != faces_to_remove.end()) { + continue; //skip removed triangles + } + Vec3i new_triangle; + for (int t_vertex_idx = 0; t_vertex_idx < 3; ++t_vertex_idx) { + size_t orig_index = mesh.indices[face_idx][t_vertex_idx]; + if (vertices_to_remove.find(orig_index) == vertices_to_remove.end()) { //this vertex was not removed + new_triangle[t_vertex_idx] = vertices_index_mapping[orig_index]; + } else { // this vertex was removed, so use the vertex mapping points to + new_triangle[t_vertex_idx] = vertices_index_mapping[vertices_index_mapping[orig_index]]; + } + } + if (new_triangle[0] == new_triangle[1] || new_triangle[1] == new_triangle[2] + || new_triangle[2] == new_triangle[0]) { + continue; //skip degenerate + } + new_indices.push_back(new_triangle); + } + + decimation_ratio = float(faces_to_remove.size()) / float(mesh.indices.size()); + +// std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; + + mesh.vertices = new_vertices; + mesh.indices = new_indices; + + vertices_index_mapping.clear(); + vertices_to_remove.clear(); + faces_to_remove.clear(); + triangles_neighbors.clear(); + face_normals.clear(); + vertex_normals.clear(); + face_indices.clear(); + triangle_count = mesh.indices.size(); + } +} + +} diff --git a/src/libslic3r/ShortEdgeCollapse.hpp b/src/libslic3r/ShortEdgeCollapse.hpp new file mode 100644 index 000000000..e6f1822c8 --- /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_ */ From 835aca60e6a5cf1b5d0b8ee22e8ff1ff6bb129fd Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Thu, 2 Jun 2022 15:04:30 +0200 Subject: [PATCH 07/11] add comments to the short edge collapse algorithm --- src/libslic3r/ShortEdgeCollapse.cpp | 81 ++++++++++++++++++----------- 1 file changed, 51 insertions(+), 30 deletions(-) diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp index 049f58b8d..0f96e93c5 100644 --- a/src/libslic3r/ShortEdgeCollapse.cpp +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -10,25 +10,31 @@ namespace Slic3r { void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) { std::unordered_map vertices_index_mapping; // this map has two usages: // in the first step, it contains mapping from original vertex index to final vertex index (which is different from original only for removed vertices) - + // in the second step, it keeps the value for the removed vertices from the first step, but for those NOT removed, it contains mapping to the new position of the new vertices vector. std::unordered_set vertices_to_remove; std::unordered_set faces_to_remove; std::vector triangles_neighbors; std::vector face_normals; - std::vector vertex_normals; - std::vector face_indices; + std::vector vertex_normals; // Vertex normal in this algorithm is normal of any! triangle that contains the vertex + + std::vector face_indices; //vector if indices, serves only the purpose of randomised traversal of the faces std::random_device rd; std::mt19937 generator(rd()); - float decimation_ratio = 1.0f; - float edge_size = 0.2f; + float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all + float edge_size = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased size_t triangle_count = mesh.indices.size(); + std::vector new_vertices; + std::vector new_indices; + + while (triangle_count > target_triangle_count) { - if (decimation_ratio < 0.4) { + // the following two switches try to keep the decimation ratio within 0.05 and 0.4 range (but nothing too clever :) + if (decimation_ratio < 0.4) { // if decimation ratio is not too large, increase it edge_size *= 1.5f; } - if (decimation_ratio < 0.05) { + if (decimation_ratio < 0.05) { // if decimation ratio is very low, increase it once more edge_size *= 1.5f; } @@ -41,81 +47,92 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ 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]; + //assign the normal of the current triangle to all its vertices. Do not care if its overridden, normal of any triangle will work + // (so in the end, each vertex has assigned normal of one of its triangles with largest face index) vertex_normals[t[0]] = n; vertex_normals[t[1]] = n; vertex_normals[t[2]] = n; - face_indices[face_idx] = face_idx; + face_indices[face_idx] = face_idx; // just an index vector, will be later shuffled for randomized face traversal } - std::vector min_vertex_score(mesh.vertices.size(), 1); + std::vector min_vertex_dot_product(mesh.vertices.size(), 1); // 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 (assigned in previous step) with the normals of faces around it. + // the lower the dot product, the more we want to keep the vertex 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_score[t[0]] = std::min(min_vertex_score[t[0]], n.dot(vertex_normals[t[0]])); - min_vertex_score[t[1]] = std::min(min_vertex_score[t[1]], n.dot(vertex_normals[t[1]])); - min_vertex_score[t[2]] = std::min(min_vertex_score[t[2]], n.dot(vertex_normals[t[2]])); + 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]])); } + // prepare vertices mapping, intitalize with self index, no vertex has been removed yet for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { vertices_index_mapping[vertex_index] = vertex_index; } + //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 (faces_to_remove.find(face_idx) != faces_to_remove.end()) { + // 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 = mesh.indices[face_idx][edge_idx]; size_t vertex_index_remove = 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; } - - if (min_vertex_score[vertex_index_remove] < min_vertex_score[vertex_index_keep]) { + // 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; } - if (min_vertex_score[vertex_index_remove] < -1.0f) { + // If vertex has already been part of a collapse, skip it (mark value -2.0 is set after collapse) + if (min_vertex_dot_product[vertex_index_remove] < -1.1f) { continue; } - + // if any of the collapsed vertices is already marked for removal, skip as well. (break, since no edge is collapsible in this case) if (vertices_to_remove.find(vertex_index_keep) != vertices_to_remove.end() || vertices_to_remove.find(vertex_index_remove) != vertices_to_remove.end()) { break; } + // find neighbour triangle over this edge. if marked for removal, skip int neighbor_face_idx = triangles_neighbors[face_idx][edge_idx]; if (neighbor_face_idx > 0 && faces_to_remove.find(neighbor_face_idx) != faces_to_remove.end()) { continue; } + //finaly do the collapse: + //mark faces for removal faces_to_remove.insert(face_idx); faces_to_remove.insert(neighbor_face_idx); + // mark one vertex for removal (with higher dot product) vertices_to_remove.insert(vertex_index_remove); + // map its index to the index of the kept vertex vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep]; - min_vertex_score[vertex_index_keep] = -2.0; + // mark the kept vertex, so that it cannot participate as a vertex for removal in any other collapse + min_vertex_dot_product[vertex_index_keep] = -2.0f; + // break, we are done with this triangle break; } } - //flatten the mapping - for (auto &pair : vertices_index_mapping) { - while (vertices_index_mapping[pair.second] != pair.second) { - pair.second = vertices_index_mapping[pair.second]; - } - } - - std::vector new_vertices; + new_vertices.clear(); new_vertices.reserve(mesh.vertices.size() - vertices_to_remove.size()); for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { + // filter out removed vertices, add those NOT removed to the new_vertices vector, and also store their new position in the index mapping if (vertices_to_remove.find(vertex_index) == vertices_to_remove.end()) { //not removed new_vertices.push_back(mesh.vertices[vertex_index]); assert(vertices_index_mapping[vertex_index] == vertex_index); @@ -123,7 +140,7 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ } } - std::vector new_indices; + new_indices.clear(); new_indices.reserve(mesh.indices.size() - faces_to_remove.size()); for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) { if (faces_to_remove.find(face_idx) != faces_to_remove.end()) { @@ -133,22 +150,26 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ for (int t_vertex_idx = 0; t_vertex_idx < 3; ++t_vertex_idx) { size_t orig_index = mesh.indices[face_idx][t_vertex_idx]; if (vertices_to_remove.find(orig_index) == vertices_to_remove.end()) { //this vertex was not removed + // NOT removed vertices point to their new position, so use directly the mapping new_triangle[t_vertex_idx] = vertices_index_mapping[orig_index]; - } else { // this vertex was removed, so use the vertex mapping points to + } else { + // Removed vertices point to the index of the kept vertex to which they collapsed. + // This vertex was surely not removed, so use the mapping twice (first to kept vertex, and then to its new position) new_triangle[t_vertex_idx] = vertices_index_mapping[vertices_index_mapping[orig_index]]; } } if (new_triangle[0] == new_triangle[1] || new_triangle[1] == new_triangle[2] || new_triangle[2] == new_triangle[0]) { - continue; //skip degenerate + continue; //skip degenerate triangles } new_indices.push_back(new_triangle); } decimation_ratio = float(faces_to_remove.size()) / float(mesh.indices.size()); -// std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; + std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; + // finally update the mesh mesh.vertices = new_vertices; mesh.indices = new_indices; From d5bf6794aacd85b75746a1eb0b722c03a1d854f6 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Thu, 2 Jun 2022 15:08:44 +0200 Subject: [PATCH 08/11] comment out debug info --- src/libslic3r/ShortEdgeCollapse.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp index 0f96e93c5..897b57759 100644 --- a/src/libslic3r/ShortEdgeCollapse.cpp +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -167,7 +167,7 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ decimation_ratio = float(faces_to_remove.size()) / float(mesh.indices.size()); - std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; +// std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; // finally update the mesh mesh.vertices = new_vertices; From 9b761d3a6f8b94e47c4c39ad60f69be771a91f77 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Thu, 2 Jun 2022 15:29:42 +0200 Subject: [PATCH 09/11] fix random generator in short edge collpase, so that results are deterministic --- src/libslic3r/ShortEdgeCollapse.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp index 897b57759..85d776843 100644 --- a/src/libslic3r/ShortEdgeCollapse.cpp +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -18,8 +18,7 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ std::vector vertex_normals; // Vertex normal in this algorithm is normal of any! triangle that contains the vertex std::vector face_indices; //vector if indices, serves only the purpose of randomised traversal of the faces - std::random_device rd; - std::mt19937 generator(rd()); + std::mt19937_64 generator { }; float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all float edge_size = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased From 13ac7a24d8429013a837c2221a74be97eb271565 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Fri, 3 Jun 2022 12:58:39 +0200 Subject: [PATCH 10/11] Refactoring of the short edge collpase, should greatly improve performance integration of NormalsUitls from SDF branch --- src/libslic3r/CMakeLists.txt | 2 + src/libslic3r/NormalUtils.cpp | 142 ++++++++++++++++ src/libslic3r/NormalUtils.hpp | 69 ++++++++ src/libslic3r/ShortEdgeCollapse.cpp | 253 ++++++++++++++-------------- 4 files changed, 338 insertions(+), 128 deletions(-) create mode 100644 src/libslic3r/NormalUtils.cpp create mode 100644 src/libslic3r/NormalUtils.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index e41d1e3cc..beabf8b02 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -166,6 +166,8 @@ add_library(libslic3r STATIC MultiPoint.cpp MultiPoint.hpp MutablePriorityQueue.hpp + NormalUtils.cpp + NormalUtils.hpp ObjectID.cpp ObjectID.hpp PerimeterGenerator.cpp diff --git a/src/libslic3r/NormalUtils.cpp b/src/libslic3r/NormalUtils.cpp new file mode 100644 index 000000000..dc9451565 --- /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 000000000..60ec57f72 --- /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/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp index 85d776843..ac55e2605 100644 --- a/src/libslic3r/ShortEdgeCollapse.cpp +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -1,4 +1,5 @@ #include "ShortEdgeCollapse.hpp" +#include "libslic3r/NormalUtils.hpp" #include #include @@ -8,57 +9,42 @@ namespace Slic3r { void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_count) { - std::unordered_map vertices_index_mapping; // this map has two usages: - // in the first step, it contains mapping from original vertex index to final vertex index (which is different from original only for removed vertices) - // in the second step, it keeps the value for the removed vertices from the first step, but for those NOT removed, it contains mapping to the new position of the new vertices vector. - std::unordered_set vertices_to_remove; - std::unordered_set faces_to_remove; - std::vector triangles_neighbors; - std::vector face_normals; - std::vector vertex_normals; // Vertex normal in this algorithm is normal of any! triangle that contains the vertex - - std::vector face_indices; //vector if indices, serves only the purpose of randomised traversal of the faces - std::mt19937_64 generator { }; - - float decimation_ratio = 1.0f; // decimation ratio updated in each iteration. it is number of removed triangles / number of all - float edge_size = 0.2f; // Allowed collapsible edge size. Starts low, but is gradually increased - size_t triangle_count = mesh.indices.size(); - - std::vector new_vertices; - std::vector new_indices; - - - while (triangle_count > target_triangle_count) { - // the following two switches try to keep the decimation ratio within 0.05 and 0.4 range (but nothing too clever :) - if (decimation_ratio < 0.4) { // if decimation ratio is not too large, increase it - edge_size *= 1.5f; + // 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]; } - if (decimation_ratio < 0.05) { // if decimation ratio is very low, increase it once more - edge_size *= 1.5f; + for (size_t i : flatten_queue) { + vertices_index_mapping[i] = idx; } + return idx; - float max_edge_len_squared = edge_size * edge_size; - triangles_neighbors = its_face_neighbors_par(mesh); - face_normals = its_face_normals(mesh); + }; - vertex_normals.resize(mesh.vertices.size()); - face_indices.resize(mesh.indices.size()); - 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]; - //assign the normal of the current triangle to all its vertices. Do not care if its overridden, normal of any triangle will work - // (so in the end, each vertex has assigned normal of one of its triangles with largest face index) - vertex_normals[t[0]] = n; - vertex_normals[t[1]] = n; - vertex_normals[t[2]] = n; + // if face is removed, mark it here + std::vector face_removal_flags(mesh.indices.size(), false); - face_indices[face_idx] = face_idx; // just an index vector, will be later shuffled for randomized face traversal - } + 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); - std::vector min_vertex_dot_product(mesh.vertices.size(), 1); // 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 (assigned in previous step) with the normals of faces around it. - // the lower the dot product, the more we want to keep the vertex 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]; @@ -66,24 +52,62 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ 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]])); } + } - // prepare vertices mapping, intitalize with self index, no vertex has been removed yet - for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { - vertices_index_mapping[vertex_index] = vertex_index; + // 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 { }; // 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 (faces_to_remove.find(face_idx) != faces_to_remove.end()) { - // if face already removed from previous collapses, skip (each collapse removes two triangles [at least] ) + 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 + + // 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 = mesh.indices[face_idx][edge_idx]; - size_t vertex_index_remove = mesh.indices[face_idx][(edge_idx + 1) % 3]; + 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) { @@ -96,91 +120,64 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ vertex_index_remove = tmp; } - // If vertex has already been part of a collapse, skip it (mark value -2.0 is set after collapse) - if (min_vertex_dot_product[vertex_index_remove] < -1.1f) { - continue; - } - // if any of the collapsed vertices is already marked for removal, skip as well. (break, since no edge is collapsible in this case) - if (vertices_to_remove.find(vertex_index_keep) != vertices_to_remove.end() - || vertices_to_remove.find(vertex_index_remove) != vertices_to_remove.end()) { - break; + //remove vertex + { + // map its index to the index of the kept vertex + vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep]; } - // find neighbour triangle over this edge. if marked for removal, skip - int neighbor_face_idx = triangles_neighbors[face_idx][edge_idx]; - if (neighbor_face_idx > 0 && faces_to_remove.find(neighbor_face_idx) != faces_to_remove.end()) { - continue; - } + 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); - //finaly do the collapse: - //mark faces for removal - faces_to_remove.insert(face_idx); - faces_to_remove.insert(neighbor_face_idx); - // mark one vertex for removal (with higher dot product) - vertices_to_remove.insert(vertex_index_remove); - // map its index to the index of the kept vertex - vertices_index_mapping[vertex_index_remove] = vertices_index_mapping[vertex_index_keep]; - // mark the kept vertex, so that it cannot participate as a vertex for removal in any other collapse - min_vertex_dot_product[vertex_index_keep] = -2.0f; - // break, we are done with this triangle + // break. this triangle is done break; } - } - new_vertices.clear(); - new_vertices.reserve(mesh.vertices.size() - vertices_to_remove.size()); - for (size_t vertex_index = 0; vertex_index < mesh.vertices.size(); ++vertex_index) { - // filter out removed vertices, add those NOT removed to the new_vertices vector, and also store their new position in the index mapping - if (vertices_to_remove.find(vertex_index) == vertices_to_remove.end()) { //not removed - new_vertices.push_back(mesh.vertices[vertex_index]); - assert(vertices_index_mapping[vertex_index] == vertex_index); - vertices_index_mapping[vertex_index] = new_vertices.size() - 1; + // 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); - new_indices.clear(); - new_indices.reserve(mesh.indices.size() - faces_to_remove.size()); - for (size_t face_idx = 0; face_idx < mesh.indices.size(); ++face_idx) { - if (faces_to_remove.find(face_idx) != faces_to_remove.end()) { - continue; //skip removed triangles - } - Vec3i new_triangle; - for (int t_vertex_idx = 0; t_vertex_idx < 3; ++t_vertex_idx) { - size_t orig_index = mesh.indices[face_idx][t_vertex_idx]; - if (vertices_to_remove.find(orig_index) == vertices_to_remove.end()) { //this vertex was not removed - // NOT removed vertices point to their new position, so use directly the mapping - new_triangle[t_vertex_idx] = vertices_index_mapping[orig_index]; - } else { - // Removed vertices point to the index of the kept vertex to which they collapsed. - // This vertex was surely not removed, so use the mapping twice (first to kept vertex, and then to its new position) - new_triangle[t_vertex_idx] = vertices_index_mapping[vertices_index_mapping[orig_index]]; - } - } - if (new_triangle[0] == new_triangle[1] || new_triangle[1] == new_triangle[2] - || new_triangle[2] == new_triangle[0]) { - continue; //skip degenerate triangles - } - new_indices.push_back(new_triangle); - } - - decimation_ratio = float(faces_to_remove.size()) / float(mesh.indices.size()); - -// std::cout << " DECIMATION RATIO: " << decimation_ratio << std::endl; - - // finally update the mesh - mesh.vertices = new_vertices; - mesh.indices = new_indices; - - vertices_index_mapping.clear(); - vertices_to_remove.clear(); - faces_to_remove.clear(); - triangles_neighbors.clear(); - face_normals.clear(); - vertex_normals.clear(); - face_indices.clear(); - triangle_count = mesh.indices.size(); + 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 + From c09781d61d809143b13cdf8567cafdb6d1bc25b4 Mon Sep 17 00:00:00 2001 From: PavelMikus Date: Mon, 6 Jun 2022 17:30:02 +0200 Subject: [PATCH 11/11] optimize embedding computation fix seed of random generators set high angle importance for nearest mode --- src/libslic3r/GCode/SeamPlacer.cpp | 23 +++++++++++++++-------- src/libslic3r/GCode/SeamPlacer.hpp | 3 ++- src/libslic3r/ShortEdgeCollapse.cpp | 2 +- src/libslic3r/TriangleSetSampling.cpp | 2 +- 4 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index 2a7448d8f..f2d386c05 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -687,8 +687,10 @@ void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) { struct SeamComparator { SeamPosition setup; + float angle_importance; 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 @@ -734,10 +736,10 @@ struct SeamComparator { // the penalites are kept close to range [0-1.x] however, it should not be relied upon float penalty_a = a.visibility + - SeamPlacer::angle_importance * compute_angle_penalty(a.local_ccw_angle) + angle_importance * compute_angle_penalty(a.local_ccw_angle) + distance_penalty_a; float penalty_b = b.visibility + - SeamPlacer::angle_importance * compute_angle_penalty(b.local_ccw_angle) + angle_importance * compute_angle_penalty(b.local_ccw_angle) + distance_penalty_b; return penalty_a < penalty_b; @@ -791,9 +793,9 @@ struct SeamComparator { } float penalty_a = a.visibility - + SeamPlacer::angle_importance * compute_angle_penalty(a.local_ccw_angle); + + angle_importance * compute_angle_penalty(a.local_ccw_angle); float penalty_b = b.visibility + - SeamPlacer::angle_importance * compute_angle_penalty(b.local_ccw_angle); + angle_importance * compute_angle_penalty(b.local_ccw_angle); return penalty_a <= penalty_b || penalty_a - penalty_b < SeamPlacer::seam_align_score_tolerance; } @@ -1071,9 +1073,14 @@ void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po) } for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { - bool layer_has_multiple_loops = - layers[layer_idx].points[0].perimeter.end_index - < layers[layer_idx].points.size() - 1; + 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; + layers[layer_idx].points[0].perimeter.end_index < layers[layer_idx].points.size() - 1; std::unique_ptr current_layer_distancer = std::make_unique(po->layers()[layer_idx]); for (SeamCandidate &perimeter_point : layers[layer_idx].points) { @@ -1082,7 +1089,7 @@ void SeamPlacer::calculate_overhangs_and_layer_embedding(const PrintObject *po) perimeter_point.overhang = prev_layer_distancer->distance_from_perimeter(point); } - if (layer_has_multiple_loops) { // search for embedded perimeter points (points hidden inside the print ,e.g. multimaterial join, best position for seam) + 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); } } diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 08974bf5f..e297497d1 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -133,7 +133,8 @@ public: static constexpr float overhang_distance_tolerance_factor = 0.5f; // determines angle importance compared to visibility ( neutral value is 1.0f. ) - static constexpr float angle_importance = 0.6f; + static constexpr float angle_importance_aligned = 0.6f; + static constexpr float angle_importance_nearest = 2.0f; // use much higher angle importance for nearest mode, to combat the visiblity info noise // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer static constexpr float enforcer_blocker_distance_tolerance = 0.35f; diff --git a/src/libslic3r/ShortEdgeCollapse.cpp b/src/libslic3r/ShortEdgeCollapse.cpp index ac55e2605..b36278c37 100644 --- a/src/libslic3r/ShortEdgeCollapse.cpp +++ b/src/libslic3r/ShortEdgeCollapse.cpp @@ -79,7 +79,7 @@ void its_short_edge_collpase(indexed_triangle_set &mesh, size_t target_triangle_ } }; - std::mt19937_64 generator { }; // default constant seed! so that results are deterministic + 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; diff --git a/src/libslic3r/TriangleSetSampling.cpp b/src/libslic3r/TriangleSetSampling.cpp index 30f9a9b40..bb03ff6d7 100644 --- a/src/libslic3r/TriangleSetSampling.cpp +++ b/src/libslic3r/TriangleSetSampling.cpp @@ -28,7 +28,7 @@ TriangleSetSamples sample_its_uniform_parallel(size_t samples_count, const index area_sum_to_triangle_idx[area_sum] = t_idx; } - std::mt19937_64 mersenne_engine { }; + std::mt19937_64 mersenne_engine { 27644437 }; // random numbers on interval [0, 1) std::uniform_real_distribution fdistribution;