diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index eb54b02dc..c97cf5010 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -11,16 +11,16 @@ #include "libslic3r/PrintConfig.hpp" #include -#include "Model.hpp" #include -#include - namespace Slic3r { namespace sla { namespace { +inline const Vec3f DOWN = {0.f, 0.f, -1.f}; +constexpr double POINTS_PER_UNIT_AREA = 1.f; + // Get the vertices of a triangle directly in an array of 3 points std::array get_triangle_vertices(const TriangleMesh &mesh, size_t faceidx) @@ -54,11 +54,11 @@ T sum_score(AccessFn &&accessfn, size_t facecount, size_t Nthreads) size_t grainsize = facecount / Nthreads; size_t from = 0, to = facecount; - return execution::reduce(ex_seq, from, to, initv, mergefn, accessfn, grainsize); + return execution::reduce(ex_tbb, from, to, initv, mergefn, accessfn, grainsize); } // Try to guess the number of support points needed to support a mesh -double get_model_supportedness(const TriangleMesh &mesh, const Transform3f &tr) +double get_misalginment_score(const TriangleMesh &mesh, const Transform3f &tr) { if (mesh.its.vertices.empty()) return std::nan(""); @@ -78,6 +78,100 @@ double get_model_supportedness(const TriangleMesh &mesh, const Transform3f &tr) return S / facecount; } +// Get area and normal of a triangle +struct Facestats { + Vec3f normal; + double area; + + explicit Facestats(const std::array &triangle) + { + Vec3f U = triangle[1] - triangle[0]; + Vec3f V = triangle[2] - triangle[0]; + Vec3f C = U.cross(V); + normal = C.normalized(); + area = 0.5 * C.norm(); + } +}; + +// The score function for a particular face +inline double get_supportedness_score(const Facestats &fc) +{ + // Simply get the angle (acos of dot product) between the face normal and + // the DOWN vector. + float phi = 1. - std::acos(fc.normal.dot(DOWN)) / float(PI); + + // Only consider faces that have have slopes below 90 deg: + phi = phi * (phi > 0.5); + + // Make the huge slopes more significant than the smaller slopes + phi = phi * phi * phi; + + // Multiply with the area of the current face + return fc.area * POINTS_PER_UNIT_AREA * phi; +} + +// Try to guess the number of support points needed to support a mesh +double get_supportedness_score(const TriangleMesh &mesh, const Transform3f &tr) +{ + if (mesh.its.vertices.empty()) return std::nan(""); + + auto accessfn = [&mesh, &tr](size_t fi) { + Facestats fc{get_transformed_triangle(mesh, tr, fi)}; + + return get_supportedness_score(fc); + }; + + size_t facecount = mesh.its.indices.size(); + size_t Nthreads = std::thread::hardware_concurrency(); + double S = unscaled(sum_score(accessfn, facecount, Nthreads)); + + return S / facecount; +} + +// Find transformed mesh ground level without copy and with parallel reduce. +float find_ground_level(const TriangleMesh &mesh, + const Transform3f & tr, + size_t threads) +{ + size_t vsize = mesh.its.vertices.size(); + + auto minfn = [](float a, float b) { return std::min(a, b); }; + + auto accessfn = [&mesh, &tr] (size_t vi) { + return (tr * mesh.its.vertices[vi]).z(); + }; + + auto zmin = std::numeric_limits::max(); + size_t granularity = vsize / threads; + return execution::reduce(ex_tbb, size_t(0), vsize, zmin, minfn, accessfn, granularity); +} + +float get_supportedness_onfloor_score(const TriangleMesh &mesh, + const Transform3f & tr) +{ + if (mesh.its.vertices.empty()) return std::nan(""); + + size_t Nthreads = std::thread::hardware_concurrency(); + + float zmin = find_ground_level(mesh, tr, Nthreads); + float zlvl = zmin + 0.1f; // Set up a slight tolerance from z level + + auto accessfn = [&mesh, &tr, zlvl](size_t fi) { + std::array tri = get_transformed_triangle(mesh, tr, fi); + Facestats fc{tri}; + + if (tri[0].z() <= zlvl && tri[1].z() <= zlvl && tri[2].z() <= zlvl) + return -fc.area * POINTS_PER_UNIT_AREA; + + return get_supportedness_score(fc); + }; + + size_t facecount = mesh.its.indices.size(); + double S = unscaled(sum_score(accessfn, facecount, Nthreads)); + + return S / facecount; +} + using XYRotation = std::array; // prepare the rotation transformation @@ -90,13 +184,107 @@ Transform3f to_transform3f(const XYRotation &rot) return rt; } +XYRotation from_transform3f(const Transform3f &tr) +{ + Vec3d rot3 = Geometry::Transformation{tr.cast()}.get_rotation(); + return {rot3.x(), rot3.y()}; +} + +inline bool is_on_floor(const SLAPrintObject &mo) +{ + auto opt_elevation = mo.config().support_object_elevation.getFloat(); + auto opt_padaround = mo.config().pad_around_object.getBool(); + + return opt_elevation < EPSILON || opt_padaround; +} + +// collect the rotations for each face of the convex hull +std::vector get_chull_rotations(const TriangleMesh &mesh, size_t max_count) +{ + TriangleMesh chull = mesh.convex_hull_3d(); + chull.require_shared_vertices(); + double chull2d_area = chull.convex_hull().area(); + double area_threshold = chull2d_area / (scaled(1e3) * scaled(1.)); + + size_t facecount = chull.its.indices.size(); + + struct RotArea { XYRotation rot; double area; }; + + auto inputs = reserve_vector(facecount); + + auto rotcmp = [](const RotArea &r1, const RotArea &r2) { + double xdiff = r1.rot[X] - r2.rot[X], ydiff = r1.rot[Y] - r2.rot[Y]; + return std::abs(xdiff) < EPSILON ? ydiff < 0. : xdiff < 0.; + }; + + auto eqcmp = [](const XYRotation &r1, const XYRotation &r2) { + double xdiff = r1[X] - r2[X], ydiff = r1[Y] - r2[Y]; + return std::abs(xdiff) < EPSILON && std::abs(ydiff) < EPSILON; + }; + + for (size_t fi = 0; fi < facecount; ++fi) { + Facestats fc{get_triangle_vertices(chull, fi)}; + + if (fc.area > area_threshold) { + auto q = Eigen::Quaternionf{}.FromTwoVectors(fc.normal, DOWN); + XYRotation rot = from_transform3f(Transform3f::Identity() * q); + RotArea ra = {rot, fc.area}; + + auto it = std::lower_bound(inputs.begin(), inputs.end(), ra, rotcmp); + + if (it == inputs.end() || !eqcmp(it->rot, rot)) + inputs.insert(it, ra); + } + } + + inputs.shrink_to_fit(); + if (!max_count) max_count = inputs.size(); + std::sort(inputs.begin(), inputs.end(), + [](const RotArea &ra, const RotArea &rb) { + return ra.area > rb.area; + }); + + auto ret = reserve_vector(std::min(max_count, inputs.size())); + for (const RotArea &ra : inputs) ret.emplace_back(ra.rot); + + return ret; +} + +// Find the best score from a set of function inputs. Evaluate for every point. +template +std::array find_min_score(Fn &&fn, It from, It to, StopCond &&stopfn) +{ + std::array ret = {}; + + double score = std::numeric_limits::max(); + + size_t Nthreads = std::thread::hardware_concurrency(); + size_t dist = std::distance(from, to); + std::vector scores(dist, score); + + execution::for_each( + ex_tbb, size_t(0), dist, [&stopfn, &scores, &fn, &from](size_t i) { + if (stopfn()) return; + + scores[i] = fn(*(from + i)); + }, + dist / Nthreads); + + auto it = std::min_element(scores.begin(), scores.end()); + + if (it != scores.end()) + ret = *(from + std::distance(scores.begin(), it)); + + return ret; +} + } // namespace -Vec2d find_best_misalignment_rotation(const SLAPrintObject & po, - float accuracy, - std::function statuscb) +Vec2d find_best_misalignment_rotation(const SLAPrintObject & po, + float accuracy, + RotOptimizeStatusCB statuscb) { - static const unsigned MAX_TRIES = 1000; + static constexpr unsigned MAX_TRIES = 1000; // return value XYRotation rot; @@ -136,22 +324,91 @@ Vec2d find_best_misalignment_rotation(const SLAPrintObject & po, // We can specify the bounds for a dimension in the following way: auto bounds = opt::bounds({ {-PI/2, PI/2}, {-PI/2, PI/2} }); - Benchmark bench; - - bench.start(); auto result = solver.to_max().optimize( [&mesh, &statusfn] (const XYRotation &rot) { statusfn(); - return get_model_supportedness(mesh, to_transform3f(rot)); + return get_misalginment_score(mesh, to_transform3f(rot)); }, opt::initvals({0., 0.}), bounds); - bench.stop(); rot = result.optimum; - std::cout << "Optimum score: " << result.score << std::endl; - std::cout << "Optimum rotation: " << result.optimum[0] << " " << result.optimum[1] << std::endl; - std::cout << "Optimization took: " << bench.getElapsedSec() << " seconds" << std::endl; + return {rot[0], rot[1]}; +} + +Vec2d find_least_supports_rotation(const SLAPrintObject & po, + float accuracy, + RotOptimizeStatusCB statuscb) +{ + static const unsigned MAX_TRIES = 1000; + + // return value + XYRotation rot; + + // We will use only one instance of this converted mesh to examine different + // rotations + TriangleMesh mesh = po.model_object()->raw_mesh(); + mesh.require_shared_vertices(); + + // To keep track of the number of iterations + unsigned status = 0; + + // The maximum number of iterations + auto max_tries = unsigned(accuracy * MAX_TRIES); + + // call status callback with zero, because we are at the start + statuscb(status); + + auto statusfn = [&statuscb, &status, &max_tries] { + // report status + statuscb(unsigned(++status * 100.0/max_tries) ); + }; + + auto stopcond = [&statuscb] { + return ! statuscb(-1); + }; + + // Different search methods have to be used depending on the model elevation + if (is_on_floor(po)) { + + std::vector inputs = get_chull_rotations(mesh, max_tries); + max_tries = inputs.size(); + + // If the model can be placed on the bed directly, we only need to + // check the 3D convex hull face rotations. + + auto objfn = [&mesh, &statusfn](const XYRotation &rot) { + statusfn(); + Transform3f tr = to_transform3f(rot); + return get_supportedness_onfloor_score(mesh, tr); + }; + + rot = find_min_score<2>(objfn, inputs.begin(), inputs.end(), stopcond); + + } else { + + // Preparing the optimizer. + size_t gridsize = std::sqrt(max_tries); // 2D grid has gridsize^2 calls + opt::Optimizer solver(opt::StopCriteria{} + .max_iterations(max_tries) + .stop_condition(stopcond), + gridsize); + + // We are searching rotations around only two axes x, y. Thus the + // problem becomes a 2 dimensional optimization task. + // We can specify the bounds for a dimension in the following way: + auto bounds = opt::bounds({ {-PI, PI}, {-PI, PI} }); + + auto result = solver.to_min().optimize( + [&mesh, &statusfn] (const XYRotation &rot) + { + statusfn(); + return get_supportedness_score(mesh, to_transform3f(rot)); + }, opt::initvals({0., 0.}), bounds); + + // Save the result + rot = result.optimum; + } return {rot[0], rot[1]}; } diff --git a/src/libslic3r/SLA/Rotfinder.hpp b/src/libslic3r/SLA/Rotfinder.hpp index 56884565f..c0007944a 100644 --- a/src/libslic3r/SLA/Rotfinder.hpp +++ b/src/libslic3r/SLA/Rotfinder.hpp @@ -37,7 +37,11 @@ Vec2d find_best_misalignment_rotation( RotOptimizeStatusCB statuscb = [] (int) { return true; } ); - +Vec2d find_least_supports_rotation( + const SLAPrintObject& modelobj, + float accuracy = 1.0f, + RotOptimizeStatusCB statuscb = [] (int) { return true; } + ); } // namespace sla } // namespace Slic3r diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp index 357036440..dfec2d6a6 100644 --- a/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp +++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp @@ -21,7 +21,7 @@ class RotoptimizeJob : public PlaterJob static inline const FindMethod Methods[] = { { L("Best misalignment"), sla::find_best_misalignment_rotation }, - { L("Least supports"), sla::find_best_misalignment_rotation } + { L("Least supports"), sla::find_least_supports_rotation } }; size_t m_method_id = 0; diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp index 59c841468..1f98463cc 100644 --- a/tests/sla_print/sla_print_tests.cpp +++ b/tests/sla_print/sla_print_tests.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include "sla_test_utils.hpp"