diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 6e0679291..6ec40f42a 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -118,10 +118,8 @@ set(SLIC3R_SOURCES GCode/PrintExtents.hpp GCode/SpiralVase.cpp GCode/SpiralVase.hpp - GCode/SeamPlacerNG.cpp - GCode/SeamPlacerNG.hpp - GCode/Subdivide.hpp - GCode/Subdivide.cpp + GCode/SeamPlacer.cpp + GCode/SeamPlacer.hpp GCode/ToolOrdering.cpp GCode/ToolOrdering.hpp GCode/WipeTower.cpp @@ -227,6 +225,8 @@ set(SLIC3R_SOURCES SlicesToTriangleMesh.cpp SlicingAdaptive.cpp SlicingAdaptive.hpp + Subdivide.cpp + Subdivide.hpp SupportMaterial.cpp SupportMaterial.hpp Surface.cpp diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp index 689897bc7..aed2ff53f 100644 --- a/src/libslic3r/GCode.cpp +++ b/src/libslic3r/GCode.cpp @@ -2312,7 +2312,6 @@ GCode::LayerResult GCode::process_layer( } // for objects // Extrude the skirt, brim, support, perimeters, infill ordered by the extruders. - std::vector> lower_layer_edge_grids(layers.size()); for (unsigned int extruder_id : layer_tools.extruders) { gcode += (layer_tools.has_wipe_tower && m_wipe_tower) ? @@ -2406,9 +2405,9 @@ GCode::LayerResult GCode::process_layer( //FIXME the following code prints regions in the order they are defined, the path is not optimized in any way. if (print.config().infill_first) { gcode += this->extrude_infill(print, by_region_specific, false); - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); } else { - gcode += this->extrude_perimeters(print, by_region_specific, lower_layer_edge_grids[instance_to_print.layer_id]); + gcode += this->extrude_perimeters(print, by_region_specific); gcode += this->extrude_infill(print,by_region_specific, false); } // ironing @@ -2543,39 +2542,11 @@ 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(); @@ -2675,14 +2646,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 ""; @@ -2703,24 +2674,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 44f4b2271..8732e797a 100644 --- a/src/libslic3r/GCode.hpp +++ b/src/libslic3r/GCode.hpp @@ -14,7 +14,7 @@ #include "GCode/SpiralVase.hpp" #include "GCode/ToolOrdering.hpp" #include "GCode/WipeTower.hpp" -#include "GCode/SeamPlacerNG.hpp" +#include "GCode/SeamPlacer.hpp" #include "GCode/GCodeProcessor.hpp" #include "EdgeGrid.hpp" #include "GCode/ThumbnailData.hpp" @@ -274,8 +274,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.); @@ -342,7 +342,7 @@ private: // For sequential print, the instance of the object to be printing has to be defined. const size_t single_object_instance_idx); - std::string extrude_perimeters(const Print &print, const std::vector &by_region, std::unique_ptr &lower_layer_edge_grid); + std::string extrude_perimeters(const Print &print, const std::vector &by_region); std::string extrude_infill(const Print &print, const std::vector &by_region, bool ironing); std::string extrude_support(const ExtrusionEntityCollection &support_fills); diff --git a/src/libslic3r/GCode/SeamPlacer.cpp b/src/libslic3r/GCode/SeamPlacer.cpp index bcd238a72..1fd7a6002 100644 --- a/src/libslic3r/GCode/SeamPlacer.cpp +++ b/src/libslic3r/GCode/SeamPlacer.cpp @@ -1,1011 +1,1327 @@ #include "SeamPlacer.hpp" +#include "tbb/parallel_for.h" +#include "tbb/blocked_range.h" +#include "tbb/parallel_reduce.h" +#include +#include +#include +#include + + +//For polynomial fitting +#include +#include + #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" + +//#define DEBUG_FILES + +#ifdef DEBUG_FILES +#include +#include +#endif namespace Slic3r { -// This penalty is added to all points inside custom blockers (subtracted from pts inside enforcers). -static constexpr float ENFORCER_BLOCKER_PENALTY = 100; +namespace SeamPlacerImpl { -// In case there are custom enforcers/blockers, the loop polygon shall always have -// sides smaller than this (so it isn't limited to original resolution). -static constexpr float MINIMAL_POLYGON_SIDE = scaled(0.2f); - -// When spAligned is active and there is a support enforcer, -// add this penalty to its center. -static constexpr float ENFORCER_CENTER_PENALTY = -10.f; - - - -// This function was introduced in 2016 to assign penalties to overhangs. -// LukasM thinks that it discriminated a bit too much, so especially external -// seams were than placed in funny places (non-overhangs were preferred too much). -// He implemented his own version (below) which applies fixed penalty for really big overlaps. -// static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) -// { -// // The extrudate is not fully supported by the lower layer. Fit a polynomial penalty curve. -// // Solved by sympy package: -// /* -// from sympy import * -// (x,a,b,c,d,r,z)=symbols('x a b c d r z') -// p = a + b*x + c*x*x + d*x*x*x -// p2 = p.subs(solve([p.subs(x, -r), p.diff(x).subs(x, -r), p.diff(x,x).subs(x, -r), p.subs(x, 0)-z], [a, b, c, d])) -// from sympy.plotting import plot -// plot(p2.subs(r,0.2).subs(z,1.), (x, -1, 3), adaptive=False, nb_of_points=400) -// */ -// if (overlap_distance < - nozzle_r) { -// // The extrudate is fully supported by the lower layer. This is the ideal case, therefore zero penalty. -// return 0.f; -// } else { -// float x = overlap_distance / nozzle_r; -// float x2 = x * x; -// float x3 = x2 * x; -// return weight_zero * (1.f + 3.f * x + 3.f * x2 + x3); -// } -// } -static float extrudate_overlap_penalty(float nozzle_r, float weight_zero, float overlap_distance) -{ - return overlap_distance > nozzle_r ? weight_zero : 0.f; +// 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); } - - -// Return a value in <0, 1> of a cubic B-spline kernel centered around zero. -// The B-spline is re-scaled so it has value 1 at zero. -static inline float bspline_kernel(float x) -{ - x = std::abs(x); - if (x < 1.f) { - return 1.f - (3.f / 2.f) * x * x + (3.f / 4.f) * x * x * x; - } - else if (x < 2.f) { - x -= 1.f; - float x2 = x * x; - float x3 = x2 * x; - return (1.f / 4.f) - (3.f / 4.f) * x + (3.f / 4.f) * x2 - (1.f / 4.f) * x3; - } - else - return 0; +Vec3f value_to_rgbf(float minimum, float maximum, float value) { + float ratio = 2.0f * (value - minimum) / (maximum - minimum); + float b = std::max(0.0f, (1.0f - ratio)); + float r = std::max(0.0f, (ratio - 1.0f)); + float g = 1.0f - b - r; + return Vec3f { r, g, b }; } - - -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(); - - Point pt_min; - double d_min = std::numeric_limits::max(); - size_t i_min = size_t(-1); - - for (size_t i = 0; i < polygon.points.size(); ++ i) { - size_t j = i + 1; - if (j == polygon.points.size()) - j = 0; - const Point &p1 = polygon.points[i]; - const Point &p2 = polygon.points[j]; - const Slic3r::Point v_seg = p2 - p1; - const Slic3r::Point v_pt = pt - p1; - const int64_t l2_seg = int64_t(v_seg(0)) * int64_t(v_seg(0)) + int64_t(v_seg(1)) * int64_t(v_seg(1)); - int64_t t_pt = int64_t(v_seg(0)) * int64_t(v_pt(0)) + int64_t(v_seg(1)) * int64_t(v_pt(1)); - if (t_pt < 0) { - // Closest to p1. - double dabs = sqrt(int64_t(v_pt(0)) * int64_t(v_pt(0)) + int64_t(v_pt(1)) * int64_t(v_pt(1))); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - pt_min = p1; - } - } - else if (t_pt > l2_seg) { - // Closest to p2. Then p2 is the starting point of another segment, which shall be discovered in the next step. - continue; - } else { - // Closest to the segment. - assert(t_pt >= 0 && t_pt <= l2_seg); - int64_t d_seg = int64_t(v_seg(1)) * int64_t(v_pt(0)) - int64_t(v_seg(0)) * int64_t(v_pt(1)); - double d = double(d_seg) / sqrt(double(l2_seg)); - double dabs = std::abs(d); - if (dabs < d_min) { - d_min = dabs; - i_min = i; - // Evaluate the foot point. - pt_min = p1; - double linv = double(d_seg) / double(l2_seg); - pt_min(0) = pt(0) - coord_t(floor(double(v_seg(1)) * linv + 0.5)); - pt_min(1) = pt(1) + coord_t(floor(double(v_seg(0)) * linv + 0.5)); - assert(Line(p1, p2).distance_to(pt_min) < scale_(1e-5)); - } - } - } - - 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; +Vec3i value_rgbi(float minimum, float maximum, float value) { + return (value_to_rgbf(minimum, maximum, value) * 255).cast(); } +//https://towardsdatascience.com/least-square-polynomial-fitting-using-c-eigen-package-c0673728bd01 +// interpolates points in z (treats z coordinates as time) and returns coefficients for axis x and y +std::vector polyfit(const std::vector &points, const std::vector &weights, size_t order) { + // check to make sure inputs are correct + assert(points.size() >= order + 1); + assert(points.size() == weights.size()); + std::vector squared_weights(weights.size()); + for (size_t index = 0; index < weights.size(); ++index) { + squared_weights[index] = sqrt(weights[index]); + } -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(); + Eigen::VectorXf V0(points.size()); + Eigen::VectorXf V1(points.size()); + Eigen::VectorXf V2(points.size()); + for (size_t index = 0; index < points.size(); index++) { + V0(index) = points[index].x() * squared_weights[index]; + V1(index) = points[index].y() * squared_weights[index]; + V2(index) = points[index].z(); + } - // 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; + // Create Matrix Placeholder of size n x k, n= number of datapoints, k = order of polynomial, for exame k = 3 for cubic polynomial + Eigen::MatrixXf T(points.size(), order + 1); + // Populate the matrix + for (size_t i = 0; i < points.size(); ++i) + { + for (size_t j = 0; j < order + 1; ++j) + { + T(i, j) = pow(V2(i), j) * squared_weights[i]; } - 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; + } + + // Solve for linear least square fit + const auto QR = T.householderQr(); + Eigen::VectorXf result0 = QR.solve(V0); + Eigen::VectorXf result1 = QR.solve(V1); + std::vector coeff { order + 1 }; + for (size_t k = 0; k < order + 1; k++) { + coeff[k] = Vec2f { result0[k], result1[k] }; + } + return coeff; +} + +Vec3f get_fitted_point(const std::vector &coefficients, float z) { + size_t order = coefficients.size() - 1; + float fitted_x = 0; + float fitted_y = 0; + for (size_t index = 0; index < order + 1; ++index) { + float z_pow = pow(z, index); + fitted_x += coefficients[index].x() * z_pow; + fitted_y += coefficients[index].y() * z_pow; + } + + return Vec3f { fitted_x, fitted_y, z }; +} + +/// Coordinate frame +class Frame { +public: + Frame() { + mX = Vec3f(1, 0, 0); + mY = Vec3f(0, 1, 0); + mZ = Vec3f(0, 0, 1); + } + + Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) : + mX(x), mY(y), mZ(z) { + } + + void set_from_z(const Vec3f &z) { + mZ = z.normalized(); + Vec3f tmpZ = mZ; + Vec3f tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3f(0, 1, 0) : Vec3f(1, 0, 0); + mY = (tmpZ.cross(tmpX)).normalized(); + mX = mY.cross(tmpZ); + } + + Vec3f to_world(const Vec3f &a) const { + return a.x() * mX + a.y() * mY + a.z() * mZ; + } + + Vec3f to_local(const Vec3f &a) const { + return Vec3f(mX.dot(a), mY.dot(a), mZ.dot(a)); + } + + const Vec3f& binormal() const { + return mX; + } + + const Vec3f& tangent() const { + return mY; + } + + const Vec3f& normal() const { + return mZ; + } + +private: + Vec3f mX, mY, mZ; +}; + +Vec3f sample_sphere_uniform(const Vec2f &samples) { + float term1 = 2.0f * float(PI) * samples.x(); + float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); + return {cos(term1) * term2, sin(term1) * term2, + 1.0f - 2.0f * samples.y()}; +} + +Vec3f sample_hemisphere_uniform(const Vec2f &samples) { + float term1 = 2.0f * float(PI) * samples.x(); + float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); + return {cos(term1) * term2, sin(term1) * term2, + abs(1.0f - 2.0f * samples.y())}; +} + +Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) { + float term1 = 2.f * float(PI) * samples.x(); + float term2 = pow(samples.y(), 1.f / (power + 1.f)); + float term3 = sqrt(1.f - term2 * term2); + + return Vec3f(cos(term1) * term3, sin(term1) * term3, term2); +} + +std::vector raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree, + const indexed_triangle_set &triangles) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start"; + + 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 }); } - while (idx_next < idx_curr && lengths.back() - lengths[idx_curr] + lengths[idx_next] < min_arm_length) - ++ idx_next; + } + + std::vector result(triangles.indices.size()); + tbb::parallel_for(tbb::blocked_range(0, result.size()), + [&](tbb::blocked_range r) { + 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); + + for (const auto &dir : precomputed_sample_directions) { + Vec3f final_ray_dir = (f.to_world(dir)); + igl::Hit hitpoint; + // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and + // direction. + Vec3d ray_origin_d = (center + normal).cast(); // start one mm above surface. + Vec3d final_ray_dir_d = final_ray_dir.cast(); + bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices, + triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint); + + if (hit) { + dest.visibility -= decrease; + } + } + } + }); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end"; + + return result; +} + +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; + } + + auto make_idx_circular = [&](int index) { + while (index < 0) { + index += polygon.size(); + } + return index % polygon.size(); + }; + + int idx_prev = 0; + int idx_curr = 0; + int 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 = make_idx_circular(idx_prev - 1); + 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 = make_idx_circular(idx_prev + 1); + } + + //push idx_next forward as far as needed + while (distance_to_next < min_arm_length) { + distance_to_next += lengths[idx_next]; + idx_next = make_idx_circular(idx_next + 1); + } + // 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; + 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(float(cross), float(dot))); + result[idx_curr] = angle; + + // increase idx_curr by one + float curr_distance = lengths[idx_curr]; + idx_curr++; + distance_to_prev += curr_distance; + distance_to_next -= curr_distance; } - return angles; + return result; } +// structure to store global information about the model - occlusion hits, enforcers, blockers +struct GlobalModelInfo { + indexed_triangle_set model; + AABBTreeIndirect::Tree<3, float> model_tree; + std::vector visiblity_info; + indexed_triangle_set enforcers; + indexed_triangle_set blockers; + AABBTreeIndirect::Tree<3, float> enforcers_tree; + AABBTreeIndirect::Tree<3, float> blockers_tree; - -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; + bool is_enforced(const Vec3f &position, float radius) const { + if (enforcers.empty()) { + return false; } - 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; - } - } - } - } - - // 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.)))); - - 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]); - - // 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; - } - } - } - } - 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; + float radius_sqr = radius * radius; + return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices, + enforcers_tree, position, radius_sqr); } - - // 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; -// } + 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); } - ++m_plan_idx; -} - - -// 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(); - 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; - } - if (! layer_po) - return last_pos; - - // Index of this layer in the respective PrintObject. - size_t layer_idx = layer_po->id() - po->layers().front()->id(); // raft layers - - assert(layer_idx < po->layer_count()); - - 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); - } - - if (seam_position != spRandom) { - // Retrieve the last start position for this object. - float last_pos_weight = 1.f; - - if (seam_position == spAligned) { - // Seam is aligned to the seam at the preceding layer. - if (po != nullptr) { - std::optional pos = m_seam_history.get_last_seam(m_po_list[po_idx], layer_idx, 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; - } - } - } - 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(); - } - 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; - break; - } - } - } - - 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 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) { + return visiblity_info[hit_idx].visibility; } 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])); + 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) { + 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) { + Points p; + perimeter->collect_points(p); + polygons.emplace_back(p); + } + } + if (polygons.empty()) { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(p); + } + } else { + Points p; + ex_entity->collect_points(p); + polygons.emplace_back(p); } } } - if (wrap_around) { - // Update first center already found. - if (out.empty()) { - // Probably an enforcer around the whole contour. Return nothing. - return out; - } - - // find last point of the enforcer at the beginning: - size_t idx = 0; - while (enforcers_idxs[idx]+1 == enforcers_idxs[idx+1]) - ++idx; - - 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); + 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 } }); } - return out; + + return polygons; } - - -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)) +// 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, std::vector &result_vec, + const GlobalModelInfo &global_model_info) { + if (orig_polygon.size() == 0) { return; - - std::vector enforcers_idxs; - std::vector blockers_idxs; - this->get_enforcers_and_blockers(layer_id, polygon, po_idx, enforcers_idxs, blockers_idxs); - - for (size_t i : enforcers_idxs) { - assert(i < penalties.size()); - penalties[i] -= float(ENFORCER_BLOCKER_PENALTY); } - for (size_t i : blockers_idxs) { - assert(i < penalties.size()); - penalties[i] += float(ENFORCER_BLOCKER_PENALTY); + + Polygon polygon = orig_polygon; + bool was_clockwise = polygon.make_counter_clockwise(); + + std::vector lengths { }; + for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) { + lengths.push_back(std::max((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm(), 0.01)); } - 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; + lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.01)); + + std::vector local_angles = calculate_polygon_angles_at_vertices(polygon, lengths, + SeamPlacer::polygon_local_angles_arm_distance); + std::shared_ptr perimeter = std::make_shared(); + + 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_vec.size(); + 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 (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { + type = EnforcedBlockedSeamPoint::Enforced; + } + + 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) + || global_model_info.is_blocked(position, distance_to_next)) { + Vec3f vec_to_next = (pos_of_next - position).normalized(); + float step_size = SeamPlacer::enforcer_blocker_oversampling_distance; + float step = step_size; + while (step < distance_to_next) { + oversampled_points.push(position + vec_to_next * step); + step += step_size; + } + } + } + + result_vec.emplace_back(position, perimeter, local_ccw_angle, type); + + } + + perimeter->end_index = result_vec.size() - 1; +} + +// 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)}; +} + +//NOTE: only rough esitmation of overhang distance +// value represents distance from edge, positive is overhang, negative is inside shape +float calculate_overhang(const SeamCandidate &point, const SeamCandidate &under_a, const SeamCandidate &under_b, + const SeamCandidate &under_c) { + auto p = Vec2d { point.position.x(), point.position.y() }; + auto a = Vec2d { under_a.position.x(), under_a.position.y() }; + auto b = Vec2d { under_b.position.x(), under_b.position.y() }; + auto c = Vec2d { under_c.position.x(), under_c.position.y() }; + + auto oriented_line_dist = [](const Vec2d a, const Vec2d b, const Vec2d p) { //signed distance from line + return -((b.x() - a.x()) * (a.y() - p.y()) - (a.x() - p.x()) * (b.y() - a.y())) / (a - b).norm(); + }; + + auto dist_ab = oriented_line_dist(a, b, p); + auto dist_bc = oriented_line_dist(b, c, p); + + // from angle and signed distances from the arms of the points on the previous layer, we + // can deduce if it is overhang and give estimation of the size. + // However, the size of the overhang is rough estimation, the sign is more reliable + if (under_b.local_ccw_angle > 0 && dist_ab > 0 && dist_bc > 0) { //convex shape, p is inside + return -((p - b).norm() + dist_ab + dist_bc) / 3.0; + } + + if (under_b.local_ccw_angle < 0 && (dist_ab < 0 || dist_bc < 0)) { //concave shape, p is inside + return -((p - b).norm() + dist_ab + dist_bc) / 3.0; + } + + return ((p - b).norm() + dist_ab + dist_bc) / 3.0; +} + +// Computes all global model info - transforms object, performs raycasting, +// stores enforces and blockers +void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB tree for raycasting and gather occlusion info: start"; +// Build AABB tree for raycasting + auto obj_transform = po->trafo_centered(); + indexed_triangle_set triangle_set; + //add all parts + for (const ModelVolume *model_volume : po->model_object()->volumes) { + if (model_volume->type() == ModelVolumeType::MODEL_PART) { + auto model_transformation = model_volume->get_matrix(); + indexed_triangle_set model_its = model_volume->mesh().its; + its_transform(model_its, model_transformation); + its_merge(triangle_set, model_its); } } -#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"); + float target_error = SeamPlacer::raycasting_decimation_target_error; + its_quadric_edge_collapse(triangle_set, 0, &target_error, nullptr, nullptr); + triangle_set = its_subdivide(triangle_set, SeamPlacer::raycasting_subdivision_target_length); + its_transform(triangle_set, obj_transform); - if (! blockers_idxs.empty()) { - ; + auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices, + triangle_set.indices); + + result.model = triangle_set; + result.model_tree = raycasting_tree; + result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB tree for raycasting and gather occlusion info: end"; + +#ifdef DEBUG_FILES + auto filename = debug_out_path(("visiblity_of_" + std::to_string(po->id().id) + ".obj").c_str()); + result.debug_export(triangle_set, filename.c_str()); +#endif +} + +void gather_enforcers_blockers(GlobalModelInfo &result, const PrintObject *po) { + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: start"; + + auto obj_transform = po->trafo(); + + for (const ModelVolume *mv : po->model_object()->volumes) { + if (mv->is_seam_painted()) { + auto model_transformation = mv->get_matrix(); + + indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER); + its_transform(enforcers, model_transformation); + its_merge(result.enforcers, enforcers); + + indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER); + its_transform(blockers, model_transformation); + its_merge(result.blockers, blockers); + } + } + its_transform(result.enforcers, obj_transform); + its_transform(result.blockers, obj_transform); + + result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices, + result.enforcers.indices); + result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices, + result.blockers.indices); + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end"; +} + +//Comparator of seam points. It has two necessary methods: is_first_better and is_first_not_much_worse +struct SeamComparator { + SeamPosition setup; + + SeamComparator(SeamPosition setup) : + setup(setup) { } + float compute_angle_penalty(float ccw_angle) const { + // This function is used: + // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x))) + // looks terribly, but it is gaussian combined with sigmoid, + // so that concave points have much smaller penalty over convex ones - size_t min_idx = std::min_element(penalties.begin(), penalties.end()) - penalties.begin(); - - for (size_t i=0; i b.type) { + return true; + } + if (b.type > a.type) { + return false; + } + + //avoid overhangs + if (a.overhang > 0.1f && b.overhang < a.overhang) { + return false; + } + + if (setup == SeamPosition::spRear) { + return a.position.y() > b.position.y(); + } + + float distance_penalty_a = 1.0f; + float distance_penalty_b = 1.0f; + if (setup == spNearest) { + distance_penalty_a = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + distance_penalty_b = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.005f); + } + + //ranges: [0 - 1] (0 - 1.3] [0.1 - 1.1) + float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle) + * distance_penalty_a; + float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(b.local_ccw_angle) + * distance_penalty_b; + + return penalty_a < penalty_b; + } + + // Comparator used during alignment. If there is close potential aligned point, it is comapred to the current + // sema point of the perimeter, to find out if the aligned point is not much worse than the current seam + bool is_first_not_much_worse(const SeamCandidate &a, const SeamCandidate &b) const { + // Blockers/Enforcers discrimination, top priority + if (a.type == EnforcedBlockedSeamPoint::Enforced) { + return true; + } + + if (a.type == EnforcedBlockedSeamPoint::Blocked) { + return false; + } + + if (a.type > b.type) { + return true; + } + if (b.type > a.type) { + return false; + } + + //avoid overhangs + if (a.overhang > 0.1f && b.overhang < a.overhang) { + return false; + } + + if (setup == SeamPosition::spRandom) { + return true; + } + + if (setup == SeamPosition::spRear) { + return a.position.y() > b.position.y(); + } + + //ranges: [0 - 1] (0 - 1.3] ; + float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle); + float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(b.local_ccw_angle); + + return penalty_a <= penalty_b || std::abs(penalty_a - penalty_b) < SeamPlacer::seam_align_score_tolerance; + } + + //always nonzero, positive + float get_penalty(const SeamCandidate &a) const { + if (setup == SeamPosition::spRear) { + return a.position.y(); + } + + return (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle); + } +} +; + +#ifdef DEBUG_FILES +void debug_export_points(const std::vector> &object_perimter_points, + const BoundingBox &bounding_box, std::string object_name, const SeamComparator &comparator) { + for (size_t layer_idx = 0; layer_idx < object_perimter_points.size(); ++layer_idx) { + std::string angles_file_name = debug_out_path((object_name + "_angles_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG angles_svg { angles_file_name, bounding_box }; + float min_vis = 0; + float max_vis = min_vis; + + float min_weight = std::numeric_limits::min(); + float max_weight = min_weight; + + for (const SeamCandidate &point : object_perimter_points[layer_idx]) { + Vec3i color = value_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.get_penalty(point)); + max_weight = std::max(max_weight, -comparator.get_penalty(point)); + + } + + std::string visiblity_file_name = debug_out_path((object_name + "_visibility_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG visibility_svg { visiblity_file_name, bounding_box }; + std::string weights_file_name = debug_out_path((object_name + "_weight_" + std::to_string(layer_idx) + ".svg").c_str()); + SVG weight_svg {weights_file_name, bounding_box }; + for (const SeamCandidate &point : object_perimter_points[layer_idx]) { + Vec3i color = value_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_rgbi(min_weight, max_weight, comparator.get_penalty(point)); + 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); + } + } +} +#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(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, he 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; + std::vector viable_indices; + std::vector viable_edges_lengths; + std::vector viable_edges; + + for (size_t index = start_index; index <= end_index; ++index) { + if (comparator.is_first_not_much_worse(perimeter_points[index], perimeter_points[viable_example_index]) && + comparator.is_first_not_much_worse(perimeter_points[viable_example_index], perimeter_points[index])) { + // index ok, push info into respective vectors + Vec3f edge_to_next; + if (index == end_index) { + edge_to_next = (perimeter_points[start_index].position - perimeter_points[index].position); + } else + { + edge_to_next = (perimeter_points[index + 1].position - perimeter_points[index].position); + } + float dist_to_next = edge_to_next.norm(); + viable_indices.push_back(index); + viable_edges_lengths.push_back(dist_to_next); + viable_edges.push_back(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; + viable_indices.clear(); + viable_edges_lengths.clear(); + viable_edges.clear(); + + Vec3f edge_to_next; + if (index == end_index) { + edge_to_next = (perimeter_points[start_index].position - perimeter_points[index].position); + } else { + edge_to_next = (perimeter_points[index + 1].position - perimeter_points[index].position); + } + float dist_to_next = edge_to_next.norm(); + viable_indices.push_back(index); + viable_edges_lengths.push_back(dist_to_next); + viable_edges.push_back(edge_to_next); + } + } + + // now pick random point from the stored options + float len_sum = std::accumulate(viable_edges_lengths.begin(), viable_edges_lengths.end(), 0.0f); + float picked_len = len_sum * (rand() / (float(RAND_MAX) + 1)); + + size_t point_idx = 0; + while (picked_len - viable_edges_lengths[point_idx] > 0) { + picked_len = picked_len - viable_edges_lengths[point_idx]; + point_idx++; + } + + Perimeter *perimeter = perimeter_points[start_index].perimeter.get(); + perimeter->seam_index = viable_indices[point_idx]; + perimeter->final_seam_position = perimeter_points[perimeter->seam_index].position + + viable_edges[point_idx].normalized() * picked_len; + perimeter->finalized = true; + +} + +} // 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 varaibles m_perimeter_points_per_object and m_perimeter_points_trees_per_object +void SeamPlacer::gather_seam_candidates(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info) { + using namespace SeamPlacerImpl; + + m_perimeter_points_per_object.emplace(po, po->layer_count()); + m_perimeter_points_trees_per_object.emplace(po, po->layer_count()); + + tbb::parallel_for(tbb::blocked_range(0, po->layers().size()), + [&](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + std::vector &layer_candidates = + m_perimeter_points_per_object[po][layer_idx]; + const Layer *layer = po->get_layer(layer_idx); + auto unscaled_z = layer->slice_z; + Polygons polygons = extract_perimeter_polygons(layer); + for (const auto &poly : polygons) { + process_perimeter_polygon(poly, unscaled_z, layer_candidates, + global_model_info); + } + auto functor = SeamCandidateCoordinateFunctor { &layer_candidates }; + m_perimeter_points_trees_per_object[po][layer_idx] = std::make_unique( + functor, layer_candidates.size()); + } + } + ); +} + +void SeamPlacer::calculate_candidates_visibility(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info) { + using namespace SeamPlacerImpl; + + tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), + [&](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + for (auto &perimeter_point : m_perimeter_points_per_object[po][layer_idx]) { + perimeter_point.visibility = global_model_info.calculate_point_visibility( + perimeter_point.position); + } + } + }); +} + +void SeamPlacer::calculate_overhangs(const PrintObject *po) { + using namespace SeamPlacerImpl; + + tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), + [&](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + for (SeamCandidate &perimeter_point : m_perimeter_points_per_object[po][layer_idx]) { + const auto calculate_layer_overhang = [&](size_t other_layer_idx) { + size_t closest_supporter = find_closest_point( + *m_perimeter_points_trees_per_object[po][other_layer_idx], + perimeter_point.position); + const SeamCandidate &supporter_point = + m_perimeter_points_per_object[po][other_layer_idx][closest_supporter]; + + auto prev_next = find_previous_and_next_perimeter_point(m_perimeter_points_per_object[po][other_layer_idx], closest_supporter); + const SeamCandidate &prev_point = + m_perimeter_points_per_object[po][other_layer_idx][prev_next.first]; + const SeamCandidate &next_point = + m_perimeter_points_per_object[po][other_layer_idx][prev_next.second]; + + return calculate_overhang(perimeter_point, prev_point, + supporter_point, next_point); + }; + + if (layer_idx > 0) { //calculate overhang + perimeter_point.overhang = calculate_layer_overhang(layer_idx-1); + } + } + } + }); + } + +// 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 +// sadly cannot be const because map access operator[] is not const, since it can create new object +bool SeamPlacer::find_next_seam_in_layer(const PrintObject *po, + std::pair &last_point_indexes, + size_t layer_idx, const SeamPlacerImpl::SeamComparator &comparator, + std::vector> &seam_string) { + using namespace SeamPlacerImpl; + + const SeamCandidate &last_point = + m_perimeter_points_per_object[po][last_point_indexes.first][last_point_indexes.second]; + + Vec3f projected_position { last_point.position.x(), last_point.position.y(), float( + po->get_layer(layer_idx)->slice_z) }; + //find closest point in next layer + size_t closest_point_index = find_closest_point( + *m_perimeter_points_trees_per_object[po][layer_idx], projected_position); + + SeamCandidate &closest_point = m_perimeter_points_per_object[po][layer_idx][closest_point_index]; + + if (closest_point.perimeter->finalized) { //already finalized, skip + return false; + } + + //from the closest point, deduce index of seam in the next layer + SeamCandidate &next_layer_seam = + m_perimeter_points_per_object[po][layer_idx][closest_point.perimeter->seam_index]; + + auto are_similar = [&](const SeamCandidate &a, const SeamCandidate &b) { + return comparator.is_first_not_much_worse(a, b) && comparator.is_first_not_much_worse(b, a); + }; + + if ((closest_point.position - projected_position).norm() < SeamPlacer::seam_align_tolerable_dist + && comparator.is_first_not_much_worse(closest_point, next_layer_seam) + && are_similar(last_point, closest_point)) { + seam_string.push_back({ layer_idx, closest_point_index }); + last_point_indexes = std::pair { layer_idx, closest_point_index }; + return true; + } else { + return false; + } + +} + +// clusters already chosen seam points into strings across multiple layers, and then +// aligns the strings via polynomial fit +// Does not change the positions of the SeamCandidates themselves, instead stores +// the new aligned position into the shared Perimeter structure of each perimeter +// Note that this position does not necesarilly lay on the perimeter. +void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) { + using namespace SeamPlacerImpl; + + // Prepares Debug files for writing. +#ifdef DEBUG_FILES + Slic3r::CNumericLocalesSetter locales_setter; + auto clusters_f = debug_out_path(("seam_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str()); + FILE *clusters = boost::nowide::fopen(clusters_f.c_str(), "w"); + if (clusters == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << clusters_f << " for writing"; + return; + } + auto aligned_f = debug_out_path(("aligned_clusters_of_" + std::to_string(po->id().id) + ".obj").c_str()); + FILE *aligns = boost::nowide::fopen(aligned_f.c_str(), "w"); + if (aligns == nullptr) { + BOOST_LOG_TRIVIAL(error) + << "stl_write_obj: Couldn't open " << clusters_f << " for writing"; + return; + } +#endif + + //gather vector of all seams on the print_object - pair of layer_index and seam__index within that layer + std::vector> seams; + for (size_t layer_idx = 0; layer_idx < m_perimeter_points_per_object[po].size(); ++layer_idx) { + std::vector &layer_perimeter_points = + m_perimeter_points_per_object[po][layer_idx]; + 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 intitializaion, this gives it better chance to choose something nice + std::sort(seams.begin(), seams.end(), + [&](const std::pair &left, const std::pair &right) { + return comparator.is_first_better(m_perimeter_points_per_object[po][left.first][left.second], + m_perimeter_points_per_object[po][right.first][right.second]); + } + ); + + //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment + for (const std::pair &seam : seams) { + size_t layer_idx = seam.first; + size_t seam_index = seam.second; + std::vector &layer_perimeter_points = + m_perimeter_points_per_object[po][layer_idx]; + if (layer_perimeter_points[seam_index].perimeter->finalized) { + // This perimeter is already aligned, skip seam + continue; + } else { + + //initialize searching for seam string - cluster of nearby seams on previous and next layers + int skips = SeamPlacer::seam_align_tolerable_skips / 2; + int next_layer = layer_idx + 1; + std::pair last_point_indexes = std::pair(layer_idx, seam_index); + + std::vector> seam_string { std::pair(layer_idx, seam_index) }; + + //find seams or potential seams in forward direction; there is a budget of skips allowed + while (skips >= 0 && next_layer < int(m_perimeter_points_per_object[po].size())) { + if (find_next_seam_in_layer(po, last_point_indexes, next_layer, comparator, seam_string)) { + //String added, last_point_pos updated, nothing to be done + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer++; + } + + //do additional check in back direction + next_layer = layer_idx - 1; + skips = SeamPlacer::seam_align_tolerable_skips / 2; + last_point_indexes = std::pair(layer_idx, seam_index); + while (skips >= 0 && next_layer >= 0) { + if (find_next_seam_in_layer(po, last_point_indexes, next_layer, comparator,seam_string)) { + //String added, last_point_pos updated, nothing to be done + } else { + // Layer skipped, reduce number of available skips + skips--; + } + next_layer--; + } + + if (seam_string.size() < seam_align_minimum_string_seams) { + //string NOT long enough to be worth aligning, skip + continue; + } + + // String is long engouh, 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) + std::vector points(seam_string.size()); + std::vector weights(seam_string.size()); + + //init min_weight by the first point + float min_weight = -comparator.get_penalty( + m_perimeter_points_per_object[po][seam_string[0].first][seam_string[0].second]); + + // In the sorted seam_string array, point which started the alignment - the best candidate + size_t best_candidate_point_index = 0; + + //gather points positions and weights, update min_weight in each step, and find the best candidate + for (size_t index = 0; index < seam_string.size(); ++index) { + points[index] = + m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second].position; + weights[index] = -comparator.get_penalty( + m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second]); + min_weight = std::min(min_weight, weights[index]); + // find the best candidate by comparing the layer indexes + if (seam_string[index].first == layer_idx) { + best_candidate_point_index = index; + } + } + + //makes all weights positive + for (float &w : weights) { + w = w - min_weight + 0.01; + } + + //NOTE: the following commented block does polynomial line fitting of the seam string. + // pre-smoothen by Laplace +// for (size_t iteration = 0; iteration < SeamPlacer::seam_align_laplace_smoothing_iterations; ++iteration) { +// std::vector new_points(seam_string.size()); +// for (int point_index = 0; point_index < points.size(); ++point_index) { +// size_t prev_idx = point_index > 0 ? point_index - 1 : point_index; +// size_t next_idx = point_index < points.size() - 1 ? point_index + 1 : point_index; +// +// new_points[point_index] = (points[prev_idx] * weights[prev_idx] +// + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / +// (weights[prev_idx] + weights[point_index] + weights[next_idx]); +// } +// points = new_points; +// } +// + // find coefficients of polynomial fit. Z coord is treated as parameter along which to fit both X and Y coords. +// std::vector coefficients = polyfit(points, weights, 4); +// +// // 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 (const auto &pair : seam_string) { +// float current_height = m_perimeter_points_per_object[po][pair.first][pair.second].position.z(); +// Vec3f seam_pos = get_fitted_point(coefficients, current_height); +// +// Perimeter *perimeter = +// m_perimeter_points_per_object[po][pair.first][pair.second].perimeter.get(); +// perimeter->final_seam_position = seam_pos; +// perimeter->finalized = true; +// } +// +// for (Vec3f &p : points) { +// p = get_fitted_point(coefficients, p.z()); +// } + + // LaPlace smoothing iterations over the gathered points. New positions from each iteration are stored in the new_points vector + // and assigned to points at the end of iteration + for (size_t iteration = 0; iteration < SeamPlacer::seam_align_laplace_smoothing_iterations; ++iteration) { + std::vector new_points(seam_string.size()); + // start from the best candidate, and smoothen down + for (int point_index = best_candidate_point_index; point_index >= 0; --point_index) { + int prev_idx = point_index > 0 ? point_index - 1 : point_index; + size_t next_idx = point_index < int(points.size()) - 1 ? point_index + 1 : point_index; + + new_points[point_index] = (points[prev_idx] * weights[prev_idx] + + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / + (weights[prev_idx] + weights[point_index] + weights[next_idx]); + } + // smoothen up the rest of the points + for (size_t point_index = best_candidate_point_index + 1; point_index < points.size(); ++point_index) { + size_t prev_idx = point_index > 0 ? point_index - 1 : point_index; + size_t next_idx = point_index < points.size() - 1 ? point_index + 1 : point_index; + + new_points[point_index] = (points[prev_idx] * weights[prev_idx] + + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / + (weights[prev_idx] + weights[point_index] + weights[next_idx]); + } + points = new_points; + } + + // Assign smoothened posiiton to each participating perimeter and set finalized flag + for (size_t index = 0; index < seam_string.size(); ++index) { + Perimeter *perimeter = + m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second].perimeter.get(); + perimeter->final_seam_position = points[index]; + 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 = m_perimeter_points_per_object[po][seam_string[i].first][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) { + Perimeter *perimeter = + m_perimeter_points_per_object[po][seam_string[i].first][seam_string[i].second].perimeter.get(); + 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_perimeter_points_trees_per_object.clear(); + m_perimeter_points_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); + 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 : start"; + calculate_overhangs(po); + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: calculate_overhangs : end"; + + BOOST_LOG_TRIVIAL(debug) + << "SeamPlacer: pick_seam_point : start"; + //pick seam point + tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), + [&](tbb::blocked_range r) { + for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { + std::vector &layer_perimeter_points = + m_perimeter_points_per_object[po][layer_idx]; + size_t current = 0; + while (current < layer_perimeter_points.size()) { + if (configured_seam_preference == spRandom) { + pick_random_seam_point(layer_perimeter_points, current); + } else { + pick_seam_point(layer_perimeter_points, current, comparator); + } + current = layer_perimeter_points[current].perimeter->end_index + 1; + } + } + }); + 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_perimeter_points_per_object[po], po->bounding_box(), std::to_string(po->id().id), + comparator); +#endif + } +} + +void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const { + using namespace SeamPlacerImpl; + const PrintObject *po = layer->object(); +//NOTE this is necessary, since layer->id() is quite unreliable + size_t layer_index = std::max(0, int(layer->id()) - int(po->slicing_parameters().raft_layers())); + double unscaled_z = layer->slice_z; + + const auto &perimeter_points_tree = *m_perimeter_points_trees_per_object.find(po)->second[layer_index]; + const auto &perimeter_points = m_perimeter_points_per_object.find(po)->second[layer_index]; + + const Point &fp = loop.first_point(); + + Vec2f unscaled_p = unscale(fp).cast(); + size_t closest_perimeter_point_index = find_closest_point(perimeter_points_tree, + Vec3f { unscaled_p.x(), unscaled_p.y(), float(unscaled_z) }); + const Perimeter *perimeter = perimeter_points[closest_perimeter_point_index].perimeter.get(); + + size_t seam_index; + if (po->config().seam_position == spNearest) { + seam_index = pick_nearest_seam_point_index(perimeter_points, perimeter->start_index, + unscale(last_pos).cast()); + } else { + seam_index = perimeter->seam_index; } - return out; + Vec3f seam_position = perimeter_points[seam_index].position; + if (perimeter->finalized) { + seam_position = perimeter->final_seam_position; + } + Point seam_point = scaled(Vec2d { seam_position.x(), seam_position.y() }); + + if (!loop.split_at_vertex(seam_point)) +// The point is not in the original loop. +// Insert it. + loop.split_at(seam_point, true); } - - -void SeamHistory::add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb) -{ - m_data_this_layer[po].push_back({pos, island_bb});; -} - - - -void SeamHistory::clear() -{ - m_layer_id = 0; - m_data_last_layer.clear(); - m_data_this_layer.clear(); -} - - -} +} // namespace Slic3r diff --git a/src/libslic3r/GCode/SeamPlacer.hpp b/src/libslic3r/GCode/SeamPlacer.hpp index 0d72939b2..2b90c1412 100644 --- a/src/libslic3r/GCode/SeamPlacer.hpp +++ b/src/libslic3r/GCode/SeamPlacer.hpp @@ -3,12 +3,15 @@ #include #include +#include +#include #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 +19,124 @@ class PrintObject; class ExtrusionLoop; class Print; class Layer; -namespace EdgeGrid { class Grid; } +namespace EdgeGrid { +class Grid; +} -class SeamHistory { -public: - SeamHistory() { clear(); } - std::optional get_last_seam(const PrintObject* po, size_t layer_id, const BoundingBox& island_bb); - void add_seam(const PrintObject* po, const Point& pos, const BoundingBox& island_bb); - void clear(); +namespace SeamPlacerImpl { -private: - struct SeamPoint { - Point m_pos; - BoundingBox m_island_bb; - }; +struct GlobalModelInfo; +struct SeamComparator; - std::map> m_data_last_layer; - std::map> m_data_this_layer; - size_t m_layer_id; +enum class EnforcedBlockedSeamPoint { + Blocked = 0, + Neutral = 1, + Enforced = 2, }; +// struct representing single perimeter loop +struct Perimeter { + size_t start_index; + size_t end_index; //inclusive! + size_t seam_index; + // 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 than further aligned +struct SeamCandidate { + SeamCandidate(const Vec3f &pos, std::shared_ptr perimeter, + float local_ccw_angle, + EnforcedBlockedSeamPoint type) : + position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), local_ccw_angle( + local_ccw_angle), type(type) { + } + const Vec3f position; + // pointer to Perimter loop of this point. It is shared across all points of the loop + const std::shared_ptr perimeter; + float visibility; + float overhang; + float local_ccw_angle; + EnforcedBlockedSeamPoint type; +}; + +struct FaceVisibilityInfo { + float visibility; +}; + +struct SeamCandidateCoordinateFunctor { + SeamCandidateCoordinateFunctor(std::vector *seam_candidates) : + seam_candidates(seam_candidates) { + } + std::vector *seam_candidates; + float operator()(size_t index, size_t dim) const { + return seam_candidates->operator[](index).position[dim]; + } +}; +} // namespace SeamPlacerImpl class SeamPlacer { public: - void init(const Print& print); + using SeamCandidatesTree = + KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>; + static constexpr float raycasting_decimation_target_error = 1.0f; + static constexpr float raycasting_subdivision_target_length = 2.0f; + //square of number of rays per triangle + static constexpr size_t sqr_rays_per_triangle = 7; - // When perimeters are printed, first call this function with the respective - // external perimeter. SeamPlacer will find a location for its seam and remember it. - // Subsequent calls to get_seam will return this position. + // arm length used during angles computation + static constexpr float polygon_local_angles_arm_distance = 0.5f; + // increases angle importance at the cost of deacreasing visibility info importance. must be > 0 + static constexpr float additional_angle_importance = 0.3f; - void plan_perimeters(const std::vector perimeters, - const Layer& layer, SeamPosition seam_position, - Point last_pos, coordf_t nozzle_dmr, const PrintObject* po, - const EdgeGrid::Grid* lower_layer_edge_grid); + // If enforcer or blocker is closer to the seam candidate than this limit, the seam candidate is set to Blocker or Enforcer + static constexpr float enforcer_blocker_distance_tolerance = 0.3f; + // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size + static constexpr float enforcer_blocker_oversampling_distance = 0.1f; - void place_seam(ExtrusionLoop& loop, const Point& last_pos, bool external_first, double nozzle_diameter, - const EdgeGrid::Grid* lower_layer_edge_grid); - + // When searching for seam clusters for alignment: + // following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer + static constexpr float seam_align_score_tolerance = 0.5f; + // seam_align_tolerable_dist - if next layer closes point is too far away, break string + static constexpr float seam_align_tolerable_dist = 1.0f; + // if the seam of the current layer is too far away, and the closest seam candidate is not very good, layer is skipped. + // this param limits the number of allowed skips + static constexpr size_t seam_align_tolerable_skips = 4; + // minimum number of seams needed in cluster to make alignemnt happen + static constexpr size_t seam_align_minimum_string_seams = 6; + // iterations of laplace smoothing + static constexpr size_t seam_align_laplace_smoothing_iterations = 20; - using TreeType = AABBTreeIndirect::Tree<2, coord_t>; - using AlignedBoxType = Eigen::AlignedBox; + //The following data structures hold all perimeter points for all PrintObject. The structure is as follows: + // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter points of the given layer + std::unordered_map>> m_perimeter_points_per_object; + // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD tree of all points of the given layer + std::unordered_map>> m_perimeter_points_trees_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); + void calculate_candidates_visibility(const PrintObject *po, + const SeamPlacerImpl::GlobalModelInfo &global_model_info); + void calculate_overhangs(const PrintObject *po); + void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator); + bool find_next_seam_in_layer(const PrintObject *po, + std::pair &last_point_indexes, + size_t layer_idx,const SeamPlacerImpl::SeamComparator &comparator, + std::vector> &seam_string); }; - -} +} // namespace Slic3r #endif // libslic3r_SeamPlacer_hpp_ diff --git a/src/libslic3r/GCode/SeamPlacerNG.cpp b/src/libslic3r/GCode/SeamPlacerNG.cpp deleted file mode 100644 index 1442c00f6..000000000 --- a/src/libslic3r/GCode/SeamPlacerNG.cpp +++ /dev/null @@ -1,1339 +0,0 @@ -#include "SeamPlacerNG.hpp" - -#include "tbb/parallel_for.h" -#include "tbb/blocked_range.h" -#include "tbb/parallel_reduce.h" -#include -#include -#include -#include - -#include "Subdivide.hpp" - -//For polynomial fitting -#include -#include - -#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" - -#define DEBUG_FILES - -#ifdef DEBUG_FILES -#include -#include -#endif - -namespace Slic3r { - -namespace SeamPlacerImpl { - -// 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); -} - -Vec3f value_to_rgbf(float minimum, float maximum, float value) { - float ratio = 2.0f * (value - minimum) / (maximum - minimum); - float b = std::max(0.0f, (1.0f - ratio)); - float r = std::max(0.0f, (ratio - 1.0f)); - float g = 1.0f - b - r; - return Vec3f { r, g, b }; -} - -Vec3i value_rgbi(float minimum, float maximum, float value) { - return (value_to_rgbf(minimum, maximum, value) * 255).cast(); -} - -//https://towardsdatascience.com/least-square-polynomial-fitting-using-c-eigen-package-c0673728bd01 -// interpolates points in z (treats z coordinates as time) and returns coefficients for axis x and y -std::vector polyfit(const std::vector &points, const std::vector &weights, size_t order) { - // check to make sure inputs are correct - assert(points.size() >= order + 1); - assert(points.size() == weights.size()); - - std::vector squared_weights(weights.size()); - for (size_t index = 0; index < weights.size(); ++index) { - squared_weights[index] = sqrt(weights[index]); - } - - Eigen::VectorXf V0(points.size()); - Eigen::VectorXf V1(points.size()); - Eigen::VectorXf V2(points.size()); - for (size_t index = 0; index < points.size(); index++) { - V0(index) = points[index].x() * squared_weights[index]; - V1(index) = points[index].y() * squared_weights[index]; - V2(index) = points[index].z(); - } - - // Create Matrix Placeholder of size n x k, n= number of datapoints, k = order of polynomial, for exame k = 3 for cubic polynomial - Eigen::MatrixXf T(points.size(), order + 1); - // Populate the matrix - for (size_t i = 0; i < points.size(); ++i) - { - for (size_t j = 0; j < order + 1; ++j) - { - T(i, j) = pow(V2(i), j) * squared_weights[i]; - } - } - - // Solve for linear least square fit - const auto QR = T.householderQr(); - Eigen::VectorXf result0 = QR.solve(V0); - Eigen::VectorXf result1 = QR.solve(V1); - std::vector coeff { order + 1 }; - for (size_t k = 0; k < order + 1; k++) { - coeff[k] = Vec2f { result0[k], result1[k] }; - } - return coeff; -} - -Vec3f get_fitted_point(const std::vector &coefficients, float z) { - size_t order = coefficients.size() - 1; - float fitted_x = 0; - float fitted_y = 0; - for (size_t index = 0; index < order + 1; ++index) { - float z_pow = pow(z, index); - fitted_x += coefficients[index].x() * z_pow; - fitted_y += coefficients[index].y() * z_pow; - } - - return Vec3f { fitted_x, fitted_y, z }; -} - -/// Coordinate frame -class Frame { -public: - Frame() { - mX = Vec3f(1, 0, 0); - mY = Vec3f(0, 1, 0); - mZ = Vec3f(0, 0, 1); - } - - Frame(const Vec3f &x, const Vec3f &y, const Vec3f &z) : - mX(x), mY(y), mZ(z) { - } - - void set_from_z(const Vec3f &z) { - mZ = z.normalized(); - Vec3f tmpZ = mZ; - Vec3f tmpX = (std::abs(tmpZ.x()) > 0.99f) ? Vec3f(0, 1, 0) : Vec3f(1, 0, 0); - mY = (tmpZ.cross(tmpX)).normalized(); - mX = mY.cross(tmpZ); - } - - Vec3f to_world(const Vec3f &a) const { - return a.x() * mX + a.y() * mY + a.z() * mZ; - } - - Vec3f to_local(const Vec3f &a) const { - return Vec3f(mX.dot(a), mY.dot(a), mZ.dot(a)); - } - - const Vec3f& binormal() const { - return mX; - } - - const Vec3f& tangent() const { - return mY; - } - - const Vec3f& normal() const { - return mZ; - } - -private: - Vec3f mX, mY, mZ; -}; - -Vec3f sample_sphere_uniform(const Vec2f &samples) { - float term1 = 2.0f * float(PI) * samples.x(); - float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); - return {cos(term1) * term2, sin(term1) * term2, - 1.0f - 2.0f * samples.y()}; -} - -Vec3f sample_hemisphere_uniform(const Vec2f &samples) { - float term1 = 2.0f * float(PI) * samples.x(); - float term2 = 2.0f * sqrt(samples.y() - samples.y() * samples.y()); - return {cos(term1) * term2, sin(term1) * term2, - abs(1.0f - 2.0f * samples.y())}; -} - -Vec3f sample_power_cosine_hemisphere(const Vec2f &samples, float power) { - float term1 = 2.f * float(PI) * samples.x(); - float term2 = pow(samples.y(), 1.f / (power + 1.f)); - float term3 = sqrt(1.f - term2 * term2); - - return Vec3f(cos(term1) * term3, sin(term1) * term3, term2); -} - -std::vector raycast_visibility(const AABBTreeIndirect::Tree<3, float> &raycasting_tree, - const indexed_triangle_set &triangles) { - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: start"; - - 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 }); - } - } - - std::vector result(triangles.indices.size()); - tbb::parallel_for(tbb::blocked_range(0, result.size()), - [&](tbb::blocked_range r) { - 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); - - for (const auto &dir : precomputed_sample_directions) { - Vec3f final_ray_dir = (f.to_world(dir)); - igl::Hit hitpoint; - // FIXME: This AABBTTreeIndirect query will not compile for float ray origin and - // direction. - Vec3d ray_origin_d = (center + normal).cast(); // start one mm above surface. - Vec3d final_ray_dir_d = final_ray_dir.cast(); - bool hit = AABBTreeIndirect::intersect_ray_first_hit(triangles.vertices, - triangles.indices, raycasting_tree, ray_origin_d, final_ray_dir_d, hitpoint); - - if (hit) { - dest.visibility -= decrease; - } - } - } - }); - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: raycast visibility for " << triangles.indices.size() << " triangles: end"; - - return result; -} - -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; - } - - auto make_idx_circular = [&](int index) { - while (index < 0) { - index += polygon.size(); - } - return index % polygon.size(); - }; - - int idx_prev = 0; - int idx_curr = 0; - int 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 = make_idx_circular(idx_prev - 1); - 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 = make_idx_circular(idx_prev + 1); - } - - //push idx_next forward as far as needed - while (distance_to_next < min_arm_length) { - distance_to_next += lengths[idx_next]; - idx_next = make_idx_circular(idx_next + 1); - } - - // 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(float(cross), float(dot))); - result[idx_curr] = angle; - - // increase idx_curr by one - float curr_distance = lengths[idx_curr]; - idx_curr++; - distance_to_prev += curr_distance; - distance_to_next -= curr_distance; - } - - return result; -} - -// structure to store global information about the model - occlusion hits, enforcers, blockers -struct GlobalModelInfo { - indexed_triangle_set model; - AABBTreeIndirect::Tree<3, float> model_tree; - std::vector visiblity_info; - indexed_triangle_set enforcers; - indexed_triangle_set blockers; - AABBTreeIndirect::Tree<3, float> enforcers_tree; - AABBTreeIndirect::Tree<3, float> blockers_tree; - - bool is_enforced(const Vec3f &position, float radius) const { - if (enforcers.empty()) { - return false; - } - float radius_sqr = radius * radius; - return AABBTreeIndirect::is_any_triangle_in_radius(enforcers.vertices, enforcers.indices, - enforcers_tree, position, radius_sqr); - } - - bool is_blocked(const Vec3f &position, float radius) const { - if (blockers.empty()) { - return false; - } - float radius_sqr = radius * radius; - return AABBTreeIndirect::is_any_triangle_in_radius(blockers.vertices, blockers.indices, - blockers_tree, position, radius_sqr); - } - - float calculate_point_visibility(const Vec3f &position) const { - size_t hit_idx; - Vec3f hit_point; - if (AABBTreeIndirect::squared_distance_to_indexed_triangle_set(model.vertices, model.indices, model_tree, - position, hit_idx, hit_point) >= 0) { - return visiblity_info[hit_idx].visibility; - } else { - return 0.0f; - } - - } - -#ifdef DEBUG_FILES - void debug_export(const indexed_triangle_set &obj_mesh, const char *file_name) const { - indexed_triangle_set divided_mesh = obj_mesh; // subdivide(obj_mesh, SeamPlacer::considered_area_radius); - 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) { - 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) { - Points p; - perimeter->collect_points(p); - polygons.emplace_back(p); - } - } - if (polygons.empty()) { - Points p; - ex_entity->collect_points(p); - polygons.emplace_back(p); - } - } else { - Points p; - ex_entity->collect_points(p); - polygons.emplace_back(p); - } - } - } - - 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 } }); - } - - 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, std::vector &result_vec, - const GlobalModelInfo &global_model_info) { - if (orig_polygon.size() == 0) { - return; - } - - Polygon polygon = orig_polygon; - bool was_clockwise = polygon.make_counter_clockwise(); - - std::vector lengths { }; - for (size_t point_idx = 0; point_idx < polygon.size() - 1; ++point_idx) { - lengths.push_back(std::max((unscale(polygon[point_idx]) - unscale(polygon[point_idx + 1])).norm(), 0.01)); - } - lengths.push_back(std::max((unscale(polygon[0]) - unscale(polygon[polygon.size() - 1])).norm(), 0.01)); - - std::vector local_angles = calculate_polygon_angles_at_vertices(polygon, lengths, - SeamPlacer::polygon_local_angles_arm_distance); - std::shared_ptr perimeter = std::make_shared(); - - 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_vec.size(); - 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 (global_model_info.is_enforced(position, SeamPlacer::enforcer_blocker_distance_tolerance)) { - type = EnforcedBlockedSeamPoint::Enforced; - } - - 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) - || global_model_info.is_blocked(position, distance_to_next)) { - Vec3f vec_to_next = (pos_of_next - position).normalized(); - float step_size = SeamPlacer::enforcer_blocker_oversampling_distance; - float step = step_size; - while (step < distance_to_next) { - oversampled_points.push(position + vec_to_next * step); - step += step_size; - } - } - } - - result_vec.emplace_back(position, perimeter, local_ccw_angle, type); - - } - - perimeter->end_index = result_vec.size() - 1; -} - -// 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)}; -} - -//NOTE: only rough esitmation of overhang distance -// value represents distance from edge, positive is overhang, negative is inside shape -float calculate_overhang(const SeamCandidate &point, const SeamCandidate &under_a, const SeamCandidate &under_b, - const SeamCandidate &under_c) { - auto p = Vec2d { point.position.x(), point.position.y() }; - auto a = Vec2d { under_a.position.x(), under_a.position.y() }; - auto b = Vec2d { under_b.position.x(), under_b.position.y() }; - auto c = Vec2d { under_c.position.x(), under_c.position.y() }; - - auto oriented_line_dist = [](const Vec2d a, const Vec2d b, const Vec2d p) { //signed distance from line - return -((b.x() - a.x()) * (a.y() - p.y()) - (a.x() - p.x()) * (b.y() - a.y())) / (a - b).norm(); - }; - - auto dist_ab = oriented_line_dist(a, b, p); - auto dist_bc = oriented_line_dist(b, c, p); - - // from angle and signed distances from the arms of the points on the previous layer, we - // can deduce if it is overhang and give estimation of the size. - // However, the size of the overhang is rough estimation, the sign is more reliable - if (under_b.local_ccw_angle > 0 && dist_ab > 0 && dist_bc > 0) { //convex shape, p is inside - return -((p - b).norm() + dist_ab + dist_bc) / 3.0; - } - - if (under_b.local_ccw_angle < 0 && (dist_ab < 0 || dist_bc < 0)) { //concave shape, p is inside - return -((p - b).norm() + dist_ab + dist_bc) / 3.0; - } - - return ((p - b).norm() + dist_ab + dist_bc) / 3.0; -} - -// Computes all global model info - transforms object, performs raycasting, -// stores enforces and blockers -void compute_global_occlusion(GlobalModelInfo &result, const PrintObject *po) { - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: build AABB tree for raycasting and gather occlusion info: start"; -// Build AABB tree for raycasting - auto obj_transform = po->trafo_centered(); - indexed_triangle_set triangle_set; - //add all parts - for (const ModelVolume *model_volume : po->model_object()->volumes) { - if (model_volume->type() == ModelVolumeType::MODEL_PART) { - auto model_transformation = model_volume->get_matrix(); - indexed_triangle_set model_its = model_volume->mesh().its; - its_transform(model_its, model_transformation); - its_merge(triangle_set, model_its); - } - } - - float target_error = SeamPlacer::raycasting_decimation_target_error; - its_quadric_edge_collapse(triangle_set, 0, &target_error, nullptr, nullptr); - triangle_set = subdivide(triangle_set, SeamPlacer::raycasting_subdivision_target_length); - its_transform(triangle_set, obj_transform); - - auto raycasting_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(triangle_set.vertices, - triangle_set.indices); - - result.model = triangle_set; - result.model_tree = raycasting_tree; - result.visiblity_info = raycast_visibility(raycasting_tree, triangle_set); - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: build AABB tree for raycasting and gather occlusion info: end"; - -#ifdef DEBUG_FILES - auto filename = "visiblity_of_" + std::to_string(po->id().id) + ".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(); - - for (const ModelVolume *mv : po->model_object()->volumes) { - if (mv->is_seam_painted()) { - auto model_transformation = mv->get_matrix(); - - indexed_triangle_set enforcers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::ENFORCER); - its_transform(enforcers, model_transformation); - its_merge(result.enforcers, enforcers); - - indexed_triangle_set blockers = mv->seam_facets.get_facets(*mv, EnforcerBlockerType::BLOCKER); - its_transform(blockers, model_transformation); - its_merge(result.blockers, blockers); - } - } - its_transform(result.enforcers, obj_transform); - its_transform(result.blockers, obj_transform); - - result.enforcers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.enforcers.vertices, - result.enforcers.indices); - result.blockers_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(result.blockers.vertices, - result.blockers.indices); - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: build AABB trees for raycasting enforcers/blockers: end"; -} - -//Comparator of seam points. It has two necessary methods: is_first_better and is_first_not_much_worse -struct SeamComparator { - SeamPosition setup; - - SeamComparator(SeamPosition setup) : - setup(setup) { - } - - float compute_angle_penalty(float ccw_angle) const { - // This function is used: - // ((ℯ^(((1)/(x^(2)*3+1)))-1)/(ℯ-1))*1+((1)/(2+ℯ^(-x))) - // looks terribly, but it is gaussian combined with sigmoid, - // so that concave points have much smaller penalty over convex ones - - return gauss(ccw_angle, 0.0f, 1.0f, 3.0f) + - 1.0f / (2 + std::exp(-ccw_angle)); // sigmoid, which heavily favourizes concave angles - } - - // 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 { - // Blockers/Enforcers discrimination, top priority - if (a.type > b.type) { - return true; - } - if (b.type > a.type) { - return false; - } - - //avoid overhangs - if (a.overhang > 0.1f && b.overhang < a.overhang) { - return false; - } - - if (setup == SeamPosition::spRear) { - return a.position.y() > b.position.y(); - } - - float distance_penalty_a = 1.0f; - float distance_penalty_b = 1.0f; - if (setup == spNearest) { - distance_penalty_a = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.01f); - distance_penalty_b = 1.1f - gauss((a.position.head<2>() - preffered_location).norm(), 0.0f, 1.0f, 0.01f); - } - - //ranges: [0 - 1] (0 - 1.3] [0.1 - 1.1) - float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle) - * distance_penalty_a; - float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(b.local_ccw_angle) - * distance_penalty_b; - - return penalty_a < penalty_b; - } - - // Comparator used during alignment. If there is close potential aligned point, it is comapred to the current - // sema point of the perimeter, to find out if the aligned point is not much worse than the current seam - bool is_first_not_much_worse(const SeamCandidate &a, const SeamCandidate &b) const { - // Blockers/Enforcers discrimination, top priority - if (a.type == EnforcedBlockedSeamPoint::Enforced) { - return true; - } - - if (a.type == EnforcedBlockedSeamPoint::Blocked) { - return false; - } - - if (a.type > b.type) { - return true; - } - if (b.type > a.type) { - return false; - } - - //avoid overhangs - if (a.overhang > 0.1f && b.overhang < a.overhang) { - return false; - } - - if (setup == SeamPosition::spRandom) { - return true; - } - - if (setup == SeamPosition::spRear) { - return a.position.y() > b.position.y(); - } - - //ranges: [0 - 1] (0 - 1.3] ; - float penalty_a = (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle); - float penalty_b = (b.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(b.local_ccw_angle); - - return penalty_a <= penalty_b || std::abs(penalty_a - penalty_b) < SeamPlacer::seam_align_score_tolerance; - } - - //always nonzero, positive - float get_penalty(const SeamCandidate &a) const { - if (setup == SeamPosition::spRear) { - return a.position.y(); - } - - return (a.visibility + SeamPlacer::additional_angle_importance) * compute_angle_penalty(a.local_ccw_angle); - } -} -; - -#ifdef DEBUG_FILES -void debug_export_points(const std::vector> &object_perimter_points, - const BoundingBox &bounding_box, std::string object_name, const SeamComparator &comparator) { - for (size_t layer_idx = 0; layer_idx < object_perimter_points.size(); ++layer_idx) { - SVG angles_svg { object_name + "_angles_" + std::to_string(layer_idx) + ".svg", 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 : object_perimter_points[layer_idx]) { - Vec3i color = value_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.get_penalty(point)); - max_weight = std::max(max_weight, -comparator.get_penalty(point)); - - } - - SVG visibility_svg { object_name + "_visibility_" + std::to_string(layer_idx) + ".svg", bounding_box }; - SVG weight_svg { object_name + "_weight_" + std::to_string(layer_idx) + ".svg", bounding_box }; - for (const SeamCandidate &point : object_perimter_points[layer_idx]) { - Vec3i color = value_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_rgbi(min_weight, max_weight, comparator.get_penalty(point)); - 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); - } - } -} -#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(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, he 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; - std::vector viable_indices; - std::vector viable_edges_lengths; - std::vector viable_edges; - - for (size_t index = start_index; index <= end_index; ++index) { - if (comparator.is_first_not_much_worse(perimeter_points[index], perimeter_points[viable_example_index]) && - comparator.is_first_not_much_worse(perimeter_points[viable_example_index], perimeter_points[index])) { - // index ok, push info into respective vectors - Vec3f edge_to_next; - if (index == end_index) { - edge_to_next = (perimeter_points[start_index].position - perimeter_points[index].position); - } else - { - edge_to_next = (perimeter_points[index + 1].position - perimeter_points[index].position); - } - float dist_to_next = edge_to_next.norm(); - viable_indices.push_back(index); - viable_edges_lengths.push_back(dist_to_next); - viable_edges.push_back(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; - viable_indices.clear(); - viable_edges_lengths.clear(); - viable_edges.clear(); - - Vec3f edge_to_next; - if (index == end_index) { - edge_to_next = (perimeter_points[start_index].position - perimeter_points[index].position); - } else { - edge_to_next = (perimeter_points[index + 1].position - perimeter_points[index].position); - } - float dist_to_next = edge_to_next.norm(); - viable_indices.push_back(index); - viable_edges_lengths.push_back(dist_to_next); - viable_edges.push_back(edge_to_next); - } - } - - // now pick random point from the stored options - float len_sum = std::accumulate(viable_edges_lengths.begin(), viable_edges_lengths.end(), 0.0f); - float picked_len = len_sum * (rand() / (float(RAND_MAX) + 1)); - - size_t point_idx = 0; - while (picked_len - viable_edges_lengths[point_idx] > 0) { - picked_len = picked_len - viable_edges_lengths[point_idx]; - point_idx++; - } - - Perimeter *perimeter = perimeter_points[start_index].perimeter.get(); - perimeter->seam_index = viable_indices[point_idx]; - perimeter->final_seam_position = perimeter_points[perimeter->seam_index].position - + viable_edges[point_idx].normalized() * picked_len; - perimeter->finalized = true; - -} - -} // 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 varaibles m_perimeter_points_per_object and m_perimeter_points_trees_per_object -void SeamPlacer::gather_seam_candidates(const PrintObject *po, - const SeamPlacerImpl::GlobalModelInfo &global_model_info) { - using namespace SeamPlacerImpl; - - m_perimeter_points_per_object.emplace(po, po->layer_count()); - m_perimeter_points_trees_per_object.emplace(po, po->layer_count()); - - tbb::parallel_for(tbb::blocked_range(0, po->layers().size()), - [&](tbb::blocked_range r) { - for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { - std::vector &layer_candidates = - m_perimeter_points_per_object[po][layer_idx]; - const Layer *layer = po->get_layer(layer_idx); - auto unscaled_z = layer->slice_z; - Polygons polygons = extract_perimeter_polygons(layer); - for (const auto &poly : polygons) { - process_perimeter_polygon(poly, unscaled_z, layer_candidates, - global_model_info); - } - auto functor = SeamCandidateCoordinateFunctor { &layer_candidates }; - m_perimeter_points_trees_per_object[po][layer_idx] = std::make_unique( - functor, layer_candidates.size()); - } - } - ); -} - -void SeamPlacer::calculate_candidates_visibility(const PrintObject *po, - const SeamPlacerImpl::GlobalModelInfo &global_model_info) { - using namespace SeamPlacerImpl; - - tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), - [&](tbb::blocked_range r) { - for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { - for (auto &perimeter_point : m_perimeter_points_per_object[po][layer_idx]) { - perimeter_point.visibility = global_model_info.calculate_point_visibility( - perimeter_point.position); - } - } - }); -} - -void SeamPlacer::calculate_overhangs(const PrintObject *po) { - using namespace SeamPlacerImpl; - - tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), - [&](tbb::blocked_range r) { - for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { - for (SeamCandidate &perimeter_point : m_perimeter_points_per_object[po][layer_idx]) { - const auto calculate_layer_overhang = [&](size_t other_layer_idx) { - size_t closest_supporter = find_closest_point( - *m_perimeter_points_trees_per_object[po][other_layer_idx], - perimeter_point.position); - const SeamCandidate &supporter_point = - m_perimeter_points_per_object[po][other_layer_idx][closest_supporter]; - - auto prev_next = find_previous_and_next_perimeter_point(m_perimeter_points_per_object[po][other_layer_idx], closest_supporter); - const SeamCandidate &prev_point = - m_perimeter_points_per_object[po][other_layer_idx][prev_next.first]; - const SeamCandidate &next_point = - m_perimeter_points_per_object[po][other_layer_idx][prev_next.second]; - - return calculate_overhang(perimeter_point, prev_point, - supporter_point, next_point); - }; - - if (layer_idx > 0) { //calculate overhang - perimeter_point.overhang = calculate_layer_overhang(layer_idx-1); - } - } - } - }); - } - -// 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 -// sadly cannot be const because map access operator[] is not const, since it can create new object -bool SeamPlacer::find_next_seam_in_layer(const PrintObject *po, - std::pair &last_point_indexes, - size_t layer_idx, const SeamPlacerImpl::SeamComparator &comparator, - std::vector> &seam_string, - std::vector> &potential_string_seams) { - using namespace SeamPlacerImpl; - - const SeamCandidate &last_point = - m_perimeter_points_per_object[po][last_point_indexes.first][last_point_indexes.second]; - - Vec3f projected_position { last_point.position.x(), last_point.position.y(), float( - po->get_layer(layer_idx)->slice_z) }; - //find closest point in next layer - size_t closest_point_index = find_closest_point( - *m_perimeter_points_trees_per_object[po][layer_idx], projected_position); - - SeamCandidate &closest_point = m_perimeter_points_per_object[po][layer_idx][closest_point_index]; - - if (closest_point.perimeter->finalized) { //already finalized, skip - return false; - } - - //from the closest point, deduce index of seam in the next layer - SeamCandidate &next_layer_seam = - m_perimeter_points_per_object[po][layer_idx][closest_point.perimeter->seam_index]; - - auto are_similar = [&](const SeamCandidate &a, const SeamCandidate &b) { - return comparator.is_first_not_much_worse(a, b) && comparator.is_first_not_much_worse(b, a); - }; - - if ((next_layer_seam.position - projected_position).norm() - < SeamPlacer::seam_align_tolerable_dist && are_similar(last_point, next_layer_seam)) { //seam point is within limits, put in the close_by_points list - seam_string.emplace_back(layer_idx, closest_point.perimeter->seam_index); - last_point_indexes = std::pair { layer_idx, closest_point.perimeter->seam_index }; - return true; - } else if ((closest_point.position - projected_position).norm() - < SeamPlacer::seam_align_tolerable_dist - && comparator.is_first_not_much_worse(closest_point, next_layer_seam) - && are_similar(last_point, closest_point)) { - //seam point is far, but if the close point is not much worse, do not count it as a skip and add it to potential_string_seams - potential_string_seams.emplace_back(layer_idx, closest_point_index); - last_point_indexes = std::pair { layer_idx, closest_point_index }; - return true; - } else { - return false; - } - -} - -// clusters already chosen seam points into strings across multiple layers, and then -// aligns the strings via polynomial fit -// Does not change the positions of the SeamCandidates themselves, instead stores -// the new aligned position into the shared Perimeter structure of each perimeter -// Note that this position does not necesarilly lay on the perimeter. -void SeamPlacer::align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator) { - using namespace SeamPlacerImpl; - - // Prepares Debug files for writing. -#ifdef DEBUG_FILES - Slic3r::CNumericLocalesSetter locales_setter; - auto clusters_f = "seam_clusters_of_" + std::to_string(po->id().id) + ".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 = "aligned_clusters_of_" + std::to_string(po->id().id) + ".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 - std::vector> seams; - for (size_t layer_idx = 0; layer_idx < m_perimeter_points_per_object[po].size(); ++layer_idx) { - std::vector &layer_perimeter_points = - m_perimeter_points_per_object[po][layer_idx]; - 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 intitializaion, this gives it better chance to choose something nice - std::sort(seams.begin(), seams.end(), - [&](const std::pair &left, const std::pair &right) { - return comparator.is_first_better(m_perimeter_points_per_object[po][left.first][left.second], - m_perimeter_points_per_object[po][right.first][right.second]); - } - ); - - //align the seam points - start with the best, and check if they are aligned, if yes, skip, else start alignment - for (const std::pair &seam : seams) { - size_t layer_idx = seam.first; - size_t seam_index = seam.second; - std::vector &layer_perimeter_points = - m_perimeter_points_per_object[po][layer_idx]; - if (layer_perimeter_points[seam_index].perimeter->finalized) { - // This perimeter is already aligned, skip seam - continue; - } else { - - //initialize searching for seam string - cluster of nearby seams on previous and next layers - int skips = SeamPlacer::seam_align_tolerable_skips / 2; - int next_layer = layer_idx + 1; - std::pair last_point_indexes = std::pair(layer_idx, seam_index); - - std::vector> seam_string { std::pair(layer_idx, seam_index) }; - std::vector> potential_string_seams; - - //find seams or potential seams in forward direction; there is a budget of skips allowed - while (skips >= 0 && next_layer < int(m_perimeter_points_per_object[po].size())) { - if (find_next_seam_in_layer(po, last_point_indexes, next_layer, comparator, seam_string, - potential_string_seams)) { - //String added, last_point_pos updated, nothing to be done - } else { - // Layer skipped, reduce number of available skips - skips--; - } - next_layer++; - } - - //do additional check in back direction - next_layer = layer_idx - 1; - skips = SeamPlacer::seam_align_tolerable_skips / 2; - last_point_indexes = std::pair(layer_idx, seam_index); - while (skips >= 0 && next_layer >= 0) { - if (find_next_seam_in_layer(po, last_point_indexes, next_layer, comparator, - seam_string, - potential_string_seams)) { - //String added, last_point_pos updated, nothing to be done - } else { - // Layer skipped, reduce number of available skips - skips--; - } - next_layer--; - } - - if (seam_string.size() + potential_string_seams.size() < seam_align_minimum_string_seams) { - //string NOT long enough to be worth aligning, skip - continue; - } - - // String is long engouh, all string seams and potential string seams gathered, now do the alignment - // first merge potential_string_seams and seam_string - seam_string.insert(seam_string.end(), potential_string_seams.begin(), potential_string_seams.end()); - //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) - std::vector points(seam_string.size()); - std::vector weights(seam_string.size()); - - //init min_weight by the first point - float min_weight = -comparator.get_penalty( - m_perimeter_points_per_object[po][seam_string[0].first][seam_string[0].second]); - - // In the sorted seam_string array, point which started the alignment - the best candidate - size_t best_candidate_point_index = 0; - - //gather points positions and weights, update min_weight in each step, and find the best candidate - for (size_t index = 0; index < seam_string.size(); ++index) { - points[index] = - m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second].position; - weights[index] = -comparator.get_penalty( - m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second]); - min_weight = std::min(min_weight, weights[index]); - // find the best candidate by comparing the layer indexes - if (seam_string[index].first == layer_idx) { - best_candidate_point_index = index; - } - } - - //makes all weights positive - for (float &w : weights) { - w = w - min_weight + 0.01; - } - - //NOTE: the following commented block does polynomial line fitting of the seam string. - // pre-smoothen by Laplace -// for (size_t iteration = 0; iteration < SeamPlacer::seam_align_laplace_smoothing_iterations; ++iteration) { -// std::vector new_points(seam_string.size()); -// for (int point_index = 0; point_index < points.size(); ++point_index) { -// size_t prev_idx = point_index > 0 ? point_index - 1 : point_index; -// size_t next_idx = point_index < points.size() - 1 ? point_index + 1 : point_index; -// -// new_points[point_index] = (points[prev_idx] * weights[prev_idx] -// + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / -// (weights[prev_idx] + weights[point_index] + weights[next_idx]); -// } -// points = new_points; -// } -// - // find coefficients of polynomial fit. Z coord is treated as parameter along which to fit both X and Y coords. -// std::vector coefficients = polyfit(points, weights, 4); -// -// // 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 (const auto &pair : seam_string) { -// float current_height = m_perimeter_points_per_object[po][pair.first][pair.second].position.z(); -// Vec3f seam_pos = get_fitted_point(coefficients, current_height); -// -// Perimeter *perimeter = -// m_perimeter_points_per_object[po][pair.first][pair.second].perimeter.get(); -// perimeter->final_seam_position = seam_pos; -// perimeter->finalized = true; -// } -// -// for (Vec3f &p : points) { -// p = get_fitted_point(coefficients, p.z()); -// } - - // LaPlace smoothing iterations over the gathered points. New positions from each iteration are stored in the new_points vector - // and assigned to points at the end of iteration - for (size_t iteration = 0; iteration < SeamPlacer::seam_align_laplace_smoothing_iterations; ++iteration) { - std::vector new_points(seam_string.size()); - // start from the best candidate, and smoothen down - for (int point_index = best_candidate_point_index; point_index >= 0; --point_index) { - int prev_idx = point_index > 0 ? point_index - 1 : point_index; - size_t next_idx = point_index < int(points.size()) - 1 ? point_index + 1 : point_index; - - new_points[point_index] = (points[prev_idx] * weights[prev_idx] - + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / - (weights[prev_idx] + weights[point_index] + weights[next_idx]); - } - // smoothen up the rest of the points - for (size_t point_index = best_candidate_point_index + 1; point_index < points.size(); ++point_index) { - size_t prev_idx = point_index > 0 ? point_index - 1 : point_index; - size_t next_idx = point_index < points.size() - 1 ? point_index + 1 : point_index; - - new_points[point_index] = (points[prev_idx] * weights[prev_idx] - + points[point_index] * weights[point_index] + points[next_idx] * weights[next_idx]) / - (weights[prev_idx] + weights[point_index] + weights[next_idx]); - } - points = new_points; - } - - // Assign smoothened posiiton to each participating perimeter and set finalized flag - for (size_t index = 0; index < seam_string.size(); ++index) { - Perimeter *perimeter = - m_perimeter_points_per_object[po][seam_string[index].first][seam_string[index].second].perimeter.get(); - perimeter->final_seam_position = points[index]; - 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 = m_perimeter_points_per_object[po][seam_string[i].first][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) { - Perimeter *perimeter = - m_perimeter_points_per_object[po][seam_string[i].first][seam_string[i].second].perimeter.get(); - 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_perimeter_points_trees_per_object.clear(); - m_perimeter_points_per_object.clear(); - - for (const PrintObject *po : print.objects()) { - - 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); - - 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); - 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); - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: calculate_candidates_visibility : end"; - } - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: calculate_overhangs : start"; - calculate_overhangs(po); - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: calculate_overhangs : end"; - - BOOST_LOG_TRIVIAL(debug) - << "SeamPlacer: pick_seam_point : start"; - //pick seam point - tbb::parallel_for(tbb::blocked_range(0, m_perimeter_points_per_object[po].size()), - [&](tbb::blocked_range r) { - for (size_t layer_idx = r.begin(); layer_idx < r.end(); ++layer_idx) { - std::vector &layer_perimeter_points = - m_perimeter_points_per_object[po][layer_idx]; - size_t current = 0; - while (current < layer_perimeter_points.size()) { - if (configured_seam_preference == spRandom) { - pick_random_seam_point(layer_perimeter_points, current); - } else { - pick_seam_point(layer_perimeter_points, current, comparator); - } - current = layer_perimeter_points[current].perimeter->end_index + 1; - } - } - }); - 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_perimeter_points_per_object[po], po->bounding_box(), std::to_string(po->id().id), - comparator); -#endif - } -} - -void SeamPlacer::place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point &last_pos) const { - using namespace SeamPlacerImpl; - const PrintObject *po = layer->object(); -//NOTE this is necessary, since layer->id() is quite unreliable - size_t layer_index = std::max(0, int(layer->id()) - int(po->slicing_parameters().raft_layers())); - double unscaled_z = layer->slice_z; - - const auto &perimeter_points_tree = *m_perimeter_points_trees_per_object.find(po)->second[layer_index]; - const auto &perimeter_points = m_perimeter_points_per_object.find(po)->second[layer_index]; - - const Point &fp = loop.first_point(); - - Vec2f unscaled_p = unscale(fp).cast(); - size_t closest_perimeter_point_index = find_closest_point(perimeter_points_tree, - Vec3f { unscaled_p.x(), unscaled_p.y(), float(unscaled_z) }); - const Perimeter *perimeter = perimeter_points[closest_perimeter_point_index].perimeter.get(); - - size_t seam_index; - if (po->config().seam_position == spNearest) { - seam_index = pick_nearest_seam_point_index(perimeter_points, perimeter->start_index, - unscale(last_pos).cast()); - } else { - seam_index = perimeter->seam_index; - } - - Vec3f seam_position = perimeter_points[seam_index].position; - if (perimeter->finalized) { - seam_position = perimeter->final_seam_position; - } - Point seam_point = scaled(Vec2d { seam_position.x(), seam_position.y() }); - - 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/SeamPlacerNG.hpp b/src/libslic3r/GCode/SeamPlacerNG.hpp deleted file mode 100644 index 3883ba004..000000000 --- a/src/libslic3r/GCode/SeamPlacerNG.hpp +++ /dev/null @@ -1,142 +0,0 @@ -#ifndef libslic3r_SeamPlacerNG_hpp_ -#define libslic3r_SeamPlacerNG_hpp_ - -#include -#include -#include -#include - -#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 { - -class PrintObject; -class ExtrusionLoop; -class Print; -class Layer; - -namespace EdgeGrid { -class Grid; -} - -namespace SeamPlacerImpl { - -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; - - // 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 than further aligned -struct SeamCandidate { - SeamCandidate(const Vec3f &pos, std::shared_ptr perimeter, - float local_ccw_angle, - EnforcedBlockedSeamPoint type) : - position(pos), perimeter(perimeter), visibility(0.0f), overhang(0.0f), local_ccw_angle( - local_ccw_angle), type(type) { - } - const Vec3f position; - // pointer to Perimter loop of this point. It is shared across all points of the loop - const std::shared_ptr perimeter; - float visibility; - float overhang; - float local_ccw_angle; - EnforcedBlockedSeamPoint type; -}; - -struct FaceVisibilityInfo { - float visibility; -}; - -struct SeamCandidateCoordinateFunctor { - SeamCandidateCoordinateFunctor(std::vector *seam_candidates) : - seam_candidates(seam_candidates) { - } - std::vector *seam_candidates; - float operator()(size_t index, size_t dim) const { - return seam_candidates->operator[](index).position[dim]; - } -}; -} // namespace SeamPlacerImpl - -class SeamPlacer { -public: - using SeamCandidatesTree = - KDTreeIndirect<3, float, SeamPlacerImpl::SeamCandidateCoordinateFunctor>; - static constexpr float raycasting_decimation_target_error = 1.0f; - 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; - - // arm length used during angles computation - static constexpr float polygon_local_angles_arm_distance = 0.5f; - static constexpr float additional_angle_importance = 0.3f; - - // 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.3f; - // For long polygon sides, if they are close to the custom seam drawings, they are oversampled with this step size - static constexpr float enforcer_blocker_oversampling_distance = 0.1f; - - // When searching for seam clusters for alignment: - // following value describes, how much worse score can point have and still be picked into seam cluster instead of original seam point on the same layer - static constexpr float seam_align_score_tolerance = 0.5f; - // seam_align_tolerable_dist - if seam is closer to the previous seam position projected to the current layer than this value, - //it belongs automaticaly to the cluster - static constexpr float seam_align_tolerable_dist = 0.5f; - // 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 alignemnt happen - static constexpr size_t seam_align_minimum_string_seams = 6; - // iterations of laplace smoothing - static constexpr size_t seam_align_laplace_smoothing_iterations = 20; - - //The following data structures hold all perimeter points for all PrintObject. The structure is as follows: - // Map of PrintObjects (PO) -> vector of layers of PO -> vector of perimeter points of the given layer - std::unordered_map>> m_perimeter_points_per_object; - // Map of PrintObjects (PO) -> vector of layers of PO -> unique_ptr to KD tree of all points of the given layer - std::unordered_map>> m_perimeter_points_trees_per_object; - - void init(const Print &print); - - void place_seam(const Layer *layer, ExtrusionLoop &loop, bool external_first, const Point& last_pos) const; - -private: - void gather_seam_candidates(const PrintObject *po, const SeamPlacerImpl::GlobalModelInfo &global_model_info); - void calculate_candidates_visibility(const PrintObject *po, - const SeamPlacerImpl::GlobalModelInfo &global_model_info); - void calculate_overhangs(const PrintObject *po); - void align_seam_points(const PrintObject *po, const SeamPlacerImpl::SeamComparator &comparator); - bool find_next_seam_in_layer(const PrintObject *po, - std::pair &last_point, - size_t layer_idx,const SeamPlacerImpl::SeamComparator &comparator, - std::vector> &seam_strings, - std::vector> &outliers); -}; - -} // namespace Slic3r - -#endif // libslic3r_SeamPlacerNG_hpp_ diff --git a/src/libslic3r/GCode/Subdivide.cpp b/src/libslic3r/Subdivide.cpp similarity index 97% rename from src/libslic3r/GCode/Subdivide.cpp rename to src/libslic3r/Subdivide.cpp index 734b1147d..31988a851 100644 --- a/src/libslic3r/GCode/Subdivide.cpp +++ b/src/libslic3r/Subdivide.cpp @@ -3,7 +3,7 @@ namespace Slic3r{ -indexed_triangle_set subdivide( +indexed_triangle_set its_subdivide( const indexed_triangle_set &its, float max_length) { // same order as key order in Edge Divides @@ -123,7 +123,7 @@ indexed_triangle_set subdivide( 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]) { + if (count_edge_vertices % 2 == 0 && key_swap == (l[i1] < l[i2])) { --index_offset; } int sign = (vs.positive_order) ? 1 : -1; @@ -161,7 +161,7 @@ indexed_triangle_set subdivide( } } - if (index_offset < count_edge_vertices-1) { + 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) { diff --git a/src/libslic3r/GCode/Subdivide.hpp b/src/libslic3r/Subdivide.hpp similarity index 63% rename from src/libslic3r/GCode/Subdivide.hpp rename to src/libslic3r/Subdivide.hpp index 3f6f96d17..f97e4b3c5 100644 --- a/src/libslic3r/GCode/Subdivide.hpp +++ b/src/libslic3r/Subdivide.hpp @@ -5,7 +5,7 @@ namespace Slic3r { -indexed_triangle_set subdivide(const indexed_triangle_set &its, float max_length); +indexed_triangle_set its_subdivide(const indexed_triangle_set &its, float max_length); }