From b4e30cc8ad08de86fee89b947c35cd401ca9021a Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Wed, 26 Aug 2020 10:25:09 +0200 Subject: [PATCH 1/9] rotation finder experiments wip --- src/libslic3r/Optimizer.hpp | 1 + src/libslic3r/SLA/Rotfinder.cpp | 155 +++++++++++++++++-------- src/libslic3r/SLA/Rotfinder.hpp | 2 +- src/slic3r/GUI/Jobs/RotoptimizeJob.cpp | 4 +- 4 files changed, 108 insertions(+), 54 deletions(-) diff --git a/src/libslic3r/Optimizer.hpp b/src/libslic3r/Optimizer.hpp index 6495ae7ff..1c94f3c1e 100644 --- a/src/libslic3r/Optimizer.hpp +++ b/src/libslic3r/Optimizer.hpp @@ -368,6 +368,7 @@ template auto score_gradient(double s, const double (&grad)[N]) using AlgNLoptGenetic = detail::NLoptAlgComb; using AlgNLoptSubplex = detail::NLoptAlg; using AlgNLoptSimplex = detail::NLoptAlg; +using AlgNLoptDIRECT = detail::NLoptAlg; // TODO: define others if needed... diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index 81ef00e6b..b4b1fae39 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -1,33 +1,113 @@ #include #include -#include +//#include +#include #include #include +#include +#include #include "Model.hpp" namespace Slic3r { namespace sla { -std::array find_best_rotation(const ModelObject& modelobj, +double area(const Vec3d &p1, const Vec3d &p2, const Vec3d &p3) { + Vec3d a = p2 - p1; + Vec3d b = p3 - p1; + Vec3d c = a.cross(b); + return 0.5 * c.norm(); +} + +using VertexFaceMap = std::vector>; + +VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) { + std::vector> vmap(mesh.its.vertices.size()); + + size_t fi = 0; + for (const Vec3i &tri : mesh.its.indices) { + for (int vi = 0; vi < tri.size(); ++vi) { + auto from = vmap[tri(vi)].begin(), to = vmap[tri(vi)].end(); + vmap[tri(vi)].insert(std::lower_bound(from, to, fi), fi); + } + } + + return vmap; +} + +// Try to guess the number of support points needed to support a mesh +double calculate_model_supportedness(const TriangleMesh & mesh, + const VertexFaceMap &vmap, + const Transform3d & tr) +{ + static const double POINTS_PER_UNIT_AREA = 1.; + static const Vec3d DOWN = {0., 0., -1.}; + + double score = 0.; + +// double zmin = mesh.bounding_box().min.z(); + +// std::vector normals(mesh.its.indices.size(), Vec3d::Zero()); + + double zmin = 0; + for (auto & v : mesh.its.vertices) + zmin = std::min(zmin, double((tr * v.cast()).z())); + + for (size_t fi = 0; fi < mesh.its.indices.size(); ++fi) { + const auto &face = mesh.its.indices[fi]; + Vec3d p1 = tr * mesh.its.vertices[face(0)].cast(); + Vec3d p2 = tr * mesh.its.vertices[face(1)].cast(); + Vec3d p3 = tr * mesh.its.vertices[face(2)].cast(); + +// auto triang = std::array {p1, p2, p3}; +// double a = area(triang.begin(), triang.end()); + double a = area(p1, p2, p3); + + double zlvl = zmin + 0.1; + if (p1.z() <= zlvl && p2.z() <= zlvl && p3.z() <= zlvl) { + score += a * POINTS_PER_UNIT_AREA; + continue; + } + + + Eigen::Vector3d U = p2 - p1; + Eigen::Vector3d V = p3 - p1; + Vec3d N = U.cross(V).normalized(); + + double phi = std::acos(N.dot(DOWN)) / PI; + + std::cout << "area: " << a << std::endl; + + score += a * POINTS_PER_UNIT_AREA * phi; +// normals[fi] = N; + } + +// for (size_t vi = 0; vi < mesh.its.vertices.size(); ++vi) { +// const std::vector &neighbors = vmap[vi]; + +// const auto &v = mesh.its.vertices[vi]; +// Vec3d vt = tr * v.cast(); +// } + + return score; +} + +std::array find_best_rotation(const ModelObject& modelobj, float accuracy, std::function statuscb, std::function stopcond) { - using libnest2d::opt::Method; - using libnest2d::opt::bound; - using libnest2d::opt::Optimizer; - using libnest2d::opt::TOptimizer; - using libnest2d::opt::StopCriteria; - - static const unsigned MAX_TRIES = 100000; + static const unsigned MAX_TRIES = 1000000; // return value - std::array rot; + std::array rot; // We will use only one instance of this converted mesh to examine different // rotations - const TriangleMesh& mesh = modelobj.raw_mesh(); + TriangleMesh mesh = modelobj.raw_mesh(); + mesh.require_shared_vertices(); +// auto vmap = create_vertex_face_map(mesh); +// simplify_mesh(mesh); // For current iteration number unsigned status = 0; @@ -44,40 +124,15 @@ std::array find_best_rotation(const ModelObject& modelobj, // the same for subsequent iterations (status goes from 0 to 100 but // iterations can be many more) auto objfunc = [&mesh, &status, &statuscb, &stopcond, max_tries] - (double rx, double ry, double rz) + (const opt::Input<2> &in) { - const TriangleMesh& m = mesh; - // prepare the rotation transformation Transform3d rt = Transform3d::Identity(); + rt.rotate(Eigen::AngleAxisd(in[1], Vec3d::UnitY())); + rt.rotate(Eigen::AngleAxisd(in[0], Vec3d::UnitX())); - rt.rotate(Eigen::AngleAxisd(rz, Vec3d::UnitZ())); - rt.rotate(Eigen::AngleAxisd(ry, Vec3d::UnitY())); - rt.rotate(Eigen::AngleAxisd(rx, Vec3d::UnitX())); - - double score = 0; - - // For all triangles we calculate the normal and sum up the dot product - // (a scalar indicating how much are two vectors aligned) with each axis - // this will result in a value that is greater if a normal is aligned - // with all axes. If the normal is aligned than the triangle itself is - // orthogonal to the axes and that is good for print quality. - - // TODO: some applications optimize for minimum z-axis cross section - // area. The current function is only an example of how to optimize. - - // Later we can add more criteria like the number of overhangs, etc... - for(size_t i = 0; i < m.stl.facet_start.size(); i++) { - Vec3d n = m.stl.facet_start[i].normal.cast(); - - // rotate the normal with the current rotation given by the solver - n = rt * n; - - // We should score against the alignment with the reference planes - score += std::abs(n.dot(Vec3d::UnitX())); - score += std::abs(n.dot(Vec3d::UnitY())); - score += std::abs(n.dot(Vec3d::UnitZ())); - } + double score = sla::calculate_model_supportedness(mesh, {}, rt); + std::cout << score << std::endl; // report status if(!stopcond()) statuscb( unsigned(++status * 100.0/max_tries) ); @@ -86,26 +141,24 @@ std::array find_best_rotation(const ModelObject& modelobj, }; // Firing up the genetic optimizer. For now it uses the nlopt library. - StopCriteria stc; - stc.max_iterations = max_tries; - stc.relative_score_difference = 1e-3; - stc.stop_condition = stopcond; // stop when stopcond returns true - TOptimizer solver(stc); + + opt::Optimizer solver(opt::StopCriteria{} + .max_iterations(max_tries) + .rel_score_diff(1e-3) + .stop_condition(stopcond)); // We are searching rotations around the three axes x, y, z. Thus the // problem becomes a 3 dimensional optimization task. // We can specify the bounds for a dimension in the following way: - auto b = bound(-PI/2, PI/2); + auto b = opt::Bound{-PI, PI}; // Now we start the optimization process with initial angles (0, 0, 0) - auto result = solver.optimize_max(objfunc, - libnest2d::opt::initvals(0.0, 0.0, 0.0), - b, b, b); + auto result = solver.to_max().optimize(objfunc, opt::initvals({0.0, 0.0}), + opt::bounds({b, b})); // Save the result and fck off rot[0] = std::get<0>(result.optimum); rot[1] = std::get<1>(result.optimum); - rot[2] = std::get<2>(result.optimum); return rot; } diff --git a/src/libslic3r/SLA/Rotfinder.hpp b/src/libslic3r/SLA/Rotfinder.hpp index 4469f9731..583703203 100644 --- a/src/libslic3r/SLA/Rotfinder.hpp +++ b/src/libslic3r/SLA/Rotfinder.hpp @@ -25,7 +25,7 @@ namespace sla { * * @return Returns the rotations around each axis (x, y, z) */ -std::array find_best_rotation( +std::array find_best_rotation( const ModelObject& modelobj, float accuracy = 1.0f, std::function statuscb = [] (unsigned) {}, diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp index c847c84b4..3fd86b13f 100644 --- a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp +++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp @@ -18,7 +18,7 @@ void RotoptimizeJob::process() auto r = sla::find_best_rotation( *o, - .005f, + 1.f, [this](unsigned s) { if (s < 100) update_status(int(s), @@ -31,7 +31,7 @@ void RotoptimizeJob::process() if (!was_canceled()) { for(ModelInstance * oi : o->instances) { - oi->set_rotation({r[X], r[Y], r[Z]}); + oi->set_rotation({r[X], r[Y], 0.}); auto trmatrix = oi->get_transformation().get_matrix(); Polygon trchull = o->convex_hull_2d(trmatrix); From c193d7c93081e0cbe6c9510436f413fe759e2b10 Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Thu, 27 Aug 2020 23:13:05 +0200 Subject: [PATCH 2/9] Brute force optimization code, buggy yet wip wip wip refactor --- src/libslic3r/CMakeLists.txt | 3 +- .../Optimize/BruteforceOptimizer.hpp | 120 ++++++++++++ .../NLoptOptimizer.hpp} | 156 +-------------- src/libslic3r/Optimize/Optimizer.hpp | 182 ++++++++++++++++++ src/libslic3r/SLA/Concurrency.hpp | 80 +++++--- src/libslic3r/SLA/Rotfinder.cpp | 111 ++++++----- src/libslic3r/SLA/SupportTreeBuildsteps.cpp | 2 +- 7 files changed, 427 insertions(+), 227 deletions(-) create mode 100644 src/libslic3r/Optimize/BruteforceOptimizer.hpp rename src/libslic3r/{Optimizer.hpp => Optimize/NLoptOptimizer.hpp} (59%) create mode 100644 src/libslic3r/Optimize/Optimizer.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 09f75c747..263920ecb 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -215,7 +215,8 @@ add_library(libslic3r STATIC SimplifyMeshImpl.hpp SimplifyMesh.cpp MarchingSquares.hpp - Optimizer.hpp + Optimize/Optimizer.hpp + Optimize/NLoptOptimizer.hpp ${OpenVDBUtils_SOURCES} SLA/Pad.hpp SLA/Pad.cpp diff --git a/src/libslic3r/Optimize/BruteforceOptimizer.hpp b/src/libslic3r/Optimize/BruteforceOptimizer.hpp new file mode 100644 index 000000000..da4472568 --- /dev/null +++ b/src/libslic3r/Optimize/BruteforceOptimizer.hpp @@ -0,0 +1,120 @@ +#ifndef BRUTEFORCEOPTIMIZER_HPP +#define BRUTEFORCEOPTIMIZER_HPP + +#include + +namespace Slic3r { namespace opt { + +namespace detail { +// Implementing a bruteforce optimizer + +template +constexpr long num_iter(const std::array &idx, size_t gridsz) +{ + long ret = 0; + for (size_t i = 0; i < N; ++i) ret += idx[i] * std::pow(gridsz, i); + return ret; +} + +struct AlgBurteForce { + bool to_min; + StopCriteria stc; + size_t gridsz; + + AlgBurteForce(const StopCriteria &cr, size_t gs): stc{cr}, gridsz{gs} {} + + template + void run(std::array &idx, + Result &result, + const Bounds &bounds, + Fn &&fn, + Cmp &&cmp) + { + if (stc.stop_condition()) return; + + if constexpr (D < 0) { + Input inp; + + auto max_iter = stc.max_iterations(); + if (max_iter && num_iter(idx, gridsz) >= max_iter) return; + + for (size_t d = 0; d < N; ++d) { + const Bound &b = bounds[d]; + double step = (b.max() - b.min()) / (gridsz - 1); + inp[d] = b.min() + idx[d] * step; + } + + auto score = fn(inp); + if (cmp(score, result.score)) { + result.score = score; + result.optimum = inp; + } + + } else { + for (size_t i = 0; i < gridsz; ++i) { + idx[D] = i; + run(idx, result, bounds, std::forward(fn), + std::forward(cmp)); + } + } + } + + template + Result optimize(Fn&& fn, + const Input &/*initvals*/, + const Bounds& bounds) + { + std::array idx = {}; + Result result; + + if (to_min) { + result.score = std::numeric_limits::max(); + run(idx, result, bounds, std::forward(fn), + std::less{}); + } + else { + result.score = std::numeric_limits::lowest(); + run(idx, result, bounds, std::forward(fn), + std::greater{}); + } + + return result; + } +}; + +} // namespace bruteforce_detail + +using AlgBruteForce = detail::AlgBurteForce; + +template<> +class Optimizer { + AlgBruteForce m_alg; + +public: + + Optimizer(const StopCriteria &cr = {}, size_t gridsz = 100) + : m_alg{cr, gridsz} + {} + + Optimizer& to_max() { m_alg.to_min = false; return *this; } + Optimizer& to_min() { m_alg.to_min = true; return *this; } + + template + Result optimize(Func&& func, + const Input &initvals, + const Bounds& bounds) + { + return m_alg.optimize(std::forward(func), initvals, bounds); + } + + Optimizer &set_criteria(const StopCriteria &cr) + { + m_alg.stc = cr; return *this; + } + + const StopCriteria &get_criteria() const { return m_alg.stc; } +}; + +}} // namespace Slic3r::opt + +#endif // BRUTEFORCEOPTIMIZER_HPP diff --git a/src/libslic3r/Optimizer.hpp b/src/libslic3r/Optimize/NLoptOptimizer.hpp similarity index 59% rename from src/libslic3r/Optimizer.hpp rename to src/libslic3r/Optimize/NLoptOptimizer.hpp index 1c94f3c1e..826b1632a 100644 --- a/src/libslic3r/Optimizer.hpp +++ b/src/libslic3r/Optimize/NLoptOptimizer.hpp @@ -12,134 +12,11 @@ #endif #include -#include -#include -#include -#include -#include -#include + +#include namespace Slic3r { namespace opt { -// A type to hold the complete result of the optimization. -template struct Result { - int resultcode; - std::array optimum; - double score; -}; - -// An interval of possible input values for optimization -class Bound { - double m_min, m_max; - -public: - Bound(double min = std::numeric_limits::min(), - double max = std::numeric_limits::max()) - : m_min(min), m_max(max) - {} - - double min() const noexcept { return m_min; } - double max() const noexcept { return m_max; } -}; - -// Helper types for optimization function input and bounds -template using Input = std::array; -template using Bounds = std::array; - -// A type for specifying the stop criteria. Setter methods can be concatenated -class StopCriteria { - - // If the absolute value difference between two scores. - double m_abs_score_diff = std::nan(""); - - // If the relative value difference between two scores. - double m_rel_score_diff = std::nan(""); - - // Stop if this value or better is found. - double m_stop_score = std::nan(""); - - // A predicate that if evaluates to true, the optimization should terminate - // and the best result found prior to termination should be returned. - std::function m_stop_condition = [] { return false; }; - - // The max allowed number of iterations. - unsigned m_max_iterations = 0; - -public: - - StopCriteria & abs_score_diff(double val) - { - m_abs_score_diff = val; return *this; - } - - double abs_score_diff() const { return m_abs_score_diff; } - - StopCriteria & rel_score_diff(double val) - { - m_rel_score_diff = val; return *this; - } - - double rel_score_diff() const { return m_rel_score_diff; } - - StopCriteria & stop_score(double val) - { - m_stop_score = val; return *this; - } - - double stop_score() const { return m_stop_score; } - - StopCriteria & max_iterations(double val) - { - m_max_iterations = val; return *this; - } - - double max_iterations() const { return m_max_iterations; } - - template StopCriteria & stop_condition(Fn &&cond) - { - m_stop_condition = cond; return *this; - } - - bool stop_condition() { return m_stop_condition(); } -}; - -// Helper class to use optimization methods involving gradient. -template struct ScoreGradient { - double score; - std::optional> gradient; - - ScoreGradient(double s, const std::array &grad) - : score{s}, gradient{grad} - {} -}; - -// Helper to be used in static_assert. -template struct always_false { enum { value = false }; }; - -// Basic interface to optimizer object -template class Optimizer { -public: - - Optimizer(const StopCriteria &) - { - static_assert (always_false::value, - "Optimizer unimplemented for given method!"); - } - - Optimizer &to_min() { return *this; } - Optimizer &to_max() { return *this; } - Optimizer &set_criteria(const StopCriteria &) { return *this; } - StopCriteria get_criteria() const { return {}; }; - - template - Result optimize(Func&& func, - const Input &initvals, - const Bounds& bounds) { return {}; } - - // optional for randomized methods: - void seed(long /*s*/) {} -}; - namespace detail { // Helper types for NLopt algorithm selection in template contexts @@ -166,19 +43,6 @@ struct IsNLoptAlg> { template using NLoptOnly = std::enable_if_t::value, T>; -// Helper to convert C style array to std::array. The copy should be optimized -// away with modern compilers. -template auto to_arr(const T *a) -{ - std::array r; - std::copy(a, a + N, std::begin(r)); - return r; -} - -template auto to_arr(const T (&a) [N]) -{ - return to_arr(static_cast(a)); -} enum class OptDir { MIN, MAX }; // Where to optimize @@ -357,24 +221,12 @@ public: void seed(long s) { m_opt.seed(s); } }; -template Bounds bounds(const Bound (&b) [N]) { return detail::to_arr(b); } -template Input initvals(const double (&a) [N]) { return detail::to_arr(a); } -template auto score_gradient(double s, const double (&grad)[N]) -{ - return ScoreGradient(s, detail::to_arr(grad)); -} - -// Predefinded NLopt algorithms that are used in the codebase +// Predefinded NLopt algorithms using AlgNLoptGenetic = detail::NLoptAlgComb; using AlgNLoptSubplex = detail::NLoptAlg; using AlgNLoptSimplex = detail::NLoptAlg; using AlgNLoptDIRECT = detail::NLoptAlg; - -// TODO: define others if needed... - -// Helper defs for pre-crafted global and local optimizers that work well. -using DefaultGlobalOptimizer = Optimizer; -using DefaultLocalOptimizer = Optimizer; +using AlgNLoptMLSL = detail::NLoptAlg; }} // namespace Slic3r::opt diff --git a/src/libslic3r/Optimize/Optimizer.hpp b/src/libslic3r/Optimize/Optimizer.hpp new file mode 100644 index 000000000..05191eba2 --- /dev/null +++ b/src/libslic3r/Optimize/Optimizer.hpp @@ -0,0 +1,182 @@ +#ifndef OPTIMIZER_HPP +#define OPTIMIZER_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace Slic3r { namespace opt { + +// A type to hold the complete result of the optimization. +template struct Result { + int resultcode; // Method dependent + std::array optimum; + double score; +}; + +// An interval of possible input values for optimization +class Bound { + double m_min, m_max; + +public: + Bound(double min = std::numeric_limits::min(), + double max = std::numeric_limits::max()) + : m_min(min), m_max(max) + {} + + double min() const noexcept { return m_min; } + double max() const noexcept { return m_max; } +}; + +// Helper types for optimization function input and bounds +template using Input = std::array; +template using Bounds = std::array; + +// A type for specifying the stop criteria. Setter methods can be concatenated +class StopCriteria { + + // If the absolute value difference between two scores. + double m_abs_score_diff = std::nan(""); + + // If the relative value difference between two scores. + double m_rel_score_diff = std::nan(""); + + // Stop if this value or better is found. + double m_stop_score = std::nan(""); + + // A predicate that if evaluates to true, the optimization should terminate + // and the best result found prior to termination should be returned. + std::function m_stop_condition = [] { return false; }; + + // The max allowed number of iterations. + unsigned m_max_iterations = 0; + +public: + + StopCriteria & abs_score_diff(double val) + { + m_abs_score_diff = val; return *this; + } + + double abs_score_diff() const { return m_abs_score_diff; } + + StopCriteria & rel_score_diff(double val) + { + m_rel_score_diff = val; return *this; + } + + double rel_score_diff() const { return m_rel_score_diff; } + + StopCriteria & stop_score(double val) + { + m_stop_score = val; return *this; + } + + double stop_score() const { return m_stop_score; } + + StopCriteria & max_iterations(double val) + { + m_max_iterations = val; return *this; + } + + double max_iterations() const { return m_max_iterations; } + + template StopCriteria & stop_condition(Fn &&cond) + { + m_stop_condition = cond; return *this; + } + + bool stop_condition() { return m_stop_condition(); } +}; + +// Helper class to use optimization methods involving gradient. +template struct ScoreGradient { + double score; + std::optional> gradient; + + ScoreGradient(double s, const std::array &grad) + : score{s}, gradient{grad} + {} +}; + +// Helper to be used in static_assert. +template struct always_false { enum { value = false }; }; + +// Basic interface to optimizer object +template class Optimizer { +public: + + Optimizer(const StopCriteria &) + { + static_assert (always_false::value, + "Optimizer unimplemented for given method!"); + } + + // Switch optimization towards function minimum + Optimizer &to_min() { return *this; } + + // Switch optimization towards function maximum + Optimizer &to_max() { return *this; } + + // Set criteria for successive optimizations + Optimizer &set_criteria(const StopCriteria &) { return *this; } + + // Get current criteria + StopCriteria get_criteria() const { return {}; }; + + // Find function minimum or maximum for Func which has has signature: + // double(const Input &input) and input with dimension N + // + // Initial starting point can be given as the second parameter. + // + // For each dimension an interval (Bound) has to be given marking the bounds + // for that dimension. + // + // initvals have to be within the specified bounds, otherwise its undefined + // behavior. + // + // Func can return a score of type double or optionally a ScoreGradient + // class to indicate the function gradient for a optimization methods that + // make use of the gradient. + template + Result optimize(Func&& /*func*/, + const Input &/*initvals*/, + const Bounds& /*bounds*/) { return {}; } + + // optional for randomized methods: + void seed(long /*s*/) {} +}; + +namespace detail { + +// Helper to convert C style array to std::array. The copy should be optimized +// away with modern compilers. +template auto to_arr(const T *a) +{ + std::array r; + std::copy(a, a + N, std::begin(r)); + return r; +} + +template auto to_arr(const T (&a) [N]) +{ + return to_arr(static_cast(a)); +} + +} // namespace detail + +// Helper functions to create bounds, initial value +template Bounds bounds(const Bound (&b) [N]) { return detail::to_arr(b); } +template Input initvals(const double (&a) [N]) { return detail::to_arr(a); } +template auto score_gradient(double s, const double (&grad)[N]) +{ + return ScoreGradient(s, detail::to_arr(grad)); +} + +}} // namespace Slic3r::opt + +#endif // OPTIMIZER_HPP diff --git a/src/libslic3r/SLA/Concurrency.hpp b/src/libslic3r/SLA/Concurrency.hpp index 93ba8c4eb..f82d6a39e 100644 --- a/src/libslic3r/SLA/Concurrency.hpp +++ b/src/libslic3r/SLA/Concurrency.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -21,28 +22,43 @@ template<> struct _ccr using SpinningMutex = tbb::spin_mutex; using BlockingMutex = tbb::mutex; + template + static IteratorOnly loop_(const tbb::blocked_range &range, Fn &&fn) + { + for (auto &el : range) fn(el); + } + + template + static IntegerOnly loop_(const tbb::blocked_range &range, Fn &&fn) + { + for (I i = range.begin(); i < range.end(); ++i) fn(i); + } + template - static IteratorOnly for_each(It from, - It to, - Fn && fn, - size_t granularity = 1) + static void for_each(It from, It to, Fn &&fn, size_t granularity = 1) { tbb::parallel_for(tbb::blocked_range{from, to, granularity}, [&fn, from](const auto &range) { - for (auto &el : range) fn(el); + loop_(range, std::forward(fn)); }); } - template - static IntegerOnly for_each(I from, - I to, - Fn && fn, - size_t granularity = 1) + template + static T reduce(I from, + I to, + const T & init, + Fn && fn, + MergeFn &&mergefn, + size_t granularity = 1) { - tbb::parallel_for(tbb::blocked_range{from, to, granularity}, - [&fn](const auto &range) { - for (I i = range.begin(); i < range.end(); ++i) fn(i); - }); + return tbb::parallel_reduce( + tbb::blocked_range{from, to, granularity}, init, + [&](const auto &range, T subinit) { + T acc = subinit; + loop_(range, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); }); + return acc; + }, + std::forward(mergefn)); } }; @@ -55,23 +71,39 @@ public: using SpinningMutex = _Mtx; using BlockingMutex = _Mtx; - template - static IteratorOnly for_each(It from, - It to, - Fn &&fn, - size_t /* ignore granularity */ = 1) + template + static IteratorOnly loop_(It from, It to, Fn &&fn) { for (auto it = from; it != to; ++it) fn(*it); } - template - static IntegerOnly for_each(I from, - I to, - Fn &&fn, - size_t /* ignore granularity */ = 1) + template + static IntegerOnly loop_(I from, I to, Fn &&fn) { for (I i = from; i < to; ++i) fn(i); } + + template + static void for_each(It from, + It to, + Fn &&fn, + size_t /* ignore granularity */ = 1) + { + loop_(from, to, std::forward(fn)); + } + + template + static IntegerOnly reduce(I from, + I to, + const T & init, + Fn && fn, + MergeFn &&mergefn, + size_t /*granularity*/ = 1) + { + T acc = init; + loop_(from, to, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); }); + return acc; + } }; using ccr = _ccr; diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index b4b1fae39..723e50eeb 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -2,23 +2,19 @@ #include //#include -#include +#include #include +#include #include #include #include #include "Model.hpp" +#include + namespace Slic3r { namespace sla { -double area(const Vec3d &p1, const Vec3d &p2, const Vec3d &p3) { - Vec3d a = p2 - p1; - Vec3d b = p3 - p1; - Vec3d c = a.cross(b); - return 0.5 * c.norm(); -} - using VertexFaceMap = std::vector>; VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) { @@ -35,61 +31,75 @@ VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) { return vmap; } +// Find transformed mesh ground level without copy and with parallell reduce. +double find_ground_level(const TriangleMesh &mesh, + const Transform3d & tr, + size_t threads) +{ + size_t vsize = mesh.its.vertices.size(); + + auto minfn = [](double a, double b) { return std::min(a, b); }; + + auto findminz = [&mesh, &tr] (size_t vi, double submin) { + Vec3d v = tr * mesh.its.vertices[vi].template cast(); + return std::min(submin, v.z()); + }; + + double zmin = mesh.its.vertices.front().z(); + + return ccr_par::reduce(size_t(0), vsize, zmin, findminz, minfn, + vsize / threads); +} + // Try to guess the number of support points needed to support a mesh double calculate_model_supportedness(const TriangleMesh & mesh, - const VertexFaceMap &vmap, +// const VertexFaceMap &vmap, const Transform3d & tr) { - static const double POINTS_PER_UNIT_AREA = 1.; - static const Vec3d DOWN = {0., 0., -1.}; + static constexpr double POINTS_PER_UNIT_AREA = 1.; - double score = 0.; + if (mesh.its.vertices.empty()) return std::nan(""); -// double zmin = mesh.bounding_box().min.z(); + size_t Nthr = std::thread::hardware_concurrency(); + size_t facesize = mesh.its.indices.size(); -// std::vector normals(mesh.its.indices.size(), Vec3d::Zero()); + double zmin = find_ground_level(mesh, tr, Nthr); - double zmin = 0; - for (auto & v : mesh.its.vertices) - zmin = std::min(zmin, double((tr * v.cast()).z())); + auto score_mergefn = [&mesh, &tr, zmin](size_t fi, double subscore) { + + static const Vec3d DOWN = {0., 0., -1.}; - for (size_t fi = 0; fi < mesh.its.indices.size(); ++fi) { const auto &face = mesh.its.indices[fi]; - Vec3d p1 = tr * mesh.its.vertices[face(0)].cast(); - Vec3d p2 = tr * mesh.its.vertices[face(1)].cast(); - Vec3d p3 = tr * mesh.its.vertices[face(2)].cast(); + Vec3d p1 = tr * mesh.its.vertices[face(0)].template cast(); + Vec3d p2 = tr * mesh.its.vertices[face(1)].template cast(); + Vec3d p3 = tr * mesh.its.vertices[face(2)].template cast(); -// auto triang = std::array {p1, p2, p3}; -// double a = area(triang.begin(), triang.end()); - double a = area(p1, p2, p3); + Vec3d U = p2 - p1; + Vec3d V = p3 - p1; + Vec3d C = U.cross(V); + Vec3d N = C.normalized(); + double area = 0.5 * C.norm(); double zlvl = zmin + 0.1; if (p1.z() <= zlvl && p2.z() <= zlvl && p3.z() <= zlvl) { - score += a * POINTS_PER_UNIT_AREA; - continue; + // score += area * POINTS_PER_UNIT_AREA; + return subscore; } + double phi = 1. - std::acos(N.dot(DOWN)) / PI; + phi = phi * (phi > 0.5); - Eigen::Vector3d U = p2 - p1; - Eigen::Vector3d V = p3 - p1; - Vec3d N = U.cross(V).normalized(); + // std::cout << "area: " << area << std::endl; - double phi = std::acos(N.dot(DOWN)) / PI; + subscore += area * POINTS_PER_UNIT_AREA * phi; - std::cout << "area: " << a << std::endl; + return subscore; + }; - score += a * POINTS_PER_UNIT_AREA * phi; -// normals[fi] = N; - } + double score = ccr_seq::reduce(size_t(0), facesize, 0., score_mergefn, + std::plus{}, facesize / Nthr); -// for (size_t vi = 0; vi < mesh.its.vertices.size(); ++vi) { -// const std::vector &neighbors = vmap[vi]; - -// const auto &v = mesh.its.vertices[vi]; -// Vec3d vt = tr * v.cast(); -// } - - return score; + return score / mesh.its.indices.size(); } std::array find_best_rotation(const ModelObject& modelobj, @@ -97,7 +107,7 @@ std::array find_best_rotation(const ModelObject& modelobj, std::function statuscb, std::function stopcond) { - static const unsigned MAX_TRIES = 1000000; + static const unsigned MAX_TRIES = 100; // return value std::array rot; @@ -126,12 +136,14 @@ std::array find_best_rotation(const ModelObject& modelobj, auto objfunc = [&mesh, &status, &statuscb, &stopcond, max_tries] (const opt::Input<2> &in) { + std::cout << "in: " << in[0] << " " << in[1] << std::endl; + // prepare the rotation transformation Transform3d rt = Transform3d::Identity(); rt.rotate(Eigen::AngleAxisd(in[1], Vec3d::UnitY())); rt.rotate(Eigen::AngleAxisd(in[0], Vec3d::UnitX())); - double score = sla::calculate_model_supportedness(mesh, {}, rt); + double score = sla::calculate_model_supportedness(mesh, rt); std::cout << score << std::endl; // report status @@ -142,10 +154,11 @@ std::array find_best_rotation(const ModelObject& modelobj, // Firing up the genetic optimizer. For now it uses the nlopt library. - opt::Optimizer solver(opt::StopCriteria{} - .max_iterations(max_tries) - .rel_score_diff(1e-3) - .stop_condition(stopcond)); + opt::Optimizer solver(opt::StopCriteria{} + .max_iterations(max_tries) + .rel_score_diff(1e-6) + .stop_condition(stopcond), + 10 /*grid size*/); // We are searching rotations around the three axes x, y, z. Thus the // problem becomes a 3 dimensional optimization task. @@ -153,7 +166,7 @@ std::array find_best_rotation(const ModelObject& modelobj, auto b = opt::Bound{-PI, PI}; // Now we start the optimization process with initial angles (0, 0, 0) - auto result = solver.to_max().optimize(objfunc, opt::initvals({0.0, 0.0}), + auto result = solver.to_min().optimize(objfunc, opt::initvals({0.0, 0.0}), opt::bounds({b, b})); // Save the result and fck off diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp index 0e7af8d50..3c39c64e6 100644 --- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp +++ b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include namespace Slic3r { From c10ff4f503208f965219d901234baa44b6d6123f Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Mon, 31 Aug 2020 18:53:44 +0200 Subject: [PATCH 3/9] fixing optimizer and concurrency::reduce --- .../Optimize/BruteforceOptimizer.hpp | 20 +++++-- src/libslic3r/SLA/Concurrency.hpp | 58 +++++++++++++----- src/libslic3r/SLA/Rotfinder.cpp | 33 +++++------ tests/libslic3r/CMakeLists.txt | 1 + tests/libslic3r/test_optimizers.cpp | 59 +++++++++++++++++++ tests/sla_print/sla_print_tests.cpp | 12 ++++ 6 files changed, 142 insertions(+), 41 deletions(-) create mode 100644 tests/libslic3r/test_optimizers.cpp diff --git a/src/libslic3r/Optimize/BruteforceOptimizer.hpp b/src/libslic3r/Optimize/BruteforceOptimizer.hpp index da4472568..960676e7b 100644 --- a/src/libslic3r/Optimize/BruteforceOptimizer.hpp +++ b/src/libslic3r/Optimize/BruteforceOptimizer.hpp @@ -1,7 +1,7 @@ #ifndef BRUTEFORCEOPTIMIZER_HPP #define BRUTEFORCEOPTIMIZER_HPP -#include +#include namespace Slic3r { namespace opt { @@ -24,19 +24,19 @@ struct AlgBurteForce { AlgBurteForce(const StopCriteria &cr, size_t gs): stc{cr}, gridsz{gs} {} template - void run(std::array &idx, + bool run(std::array &idx, Result &result, const Bounds &bounds, Fn &&fn, Cmp &&cmp) { - if (stc.stop_condition()) return; + if (stc.stop_condition()) return false; if constexpr (D < 0) { Input inp; auto max_iter = stc.max_iterations(); - if (max_iter && num_iter(idx, gridsz) >= max_iter) return; + if (max_iter && num_iter(idx, gridsz) >= max_iter) return false; for (size_t d = 0; d < N; ++d) { const Bound &b = bounds[d]; @@ -46,17 +46,25 @@ struct AlgBurteForce { auto score = fn(inp); if (cmp(score, result.score)) { + double absdiff = std::abs(score - result.score); + result.score = score; result.optimum = inp; + + if (absdiff < stc.abs_score_diff() || + absdiff < stc.rel_score_diff() * std::abs(score)) + return false; } } else { for (size_t i = 0; i < gridsz; ++i) { idx[D] = i; - run(idx, result, bounds, std::forward(fn), - std::forward(cmp)); + if (!run(idx, result, bounds, std::forward(fn), + std::forward(cmp))) return false; } } + + return true; } template diff --git a/src/libslic3r/SLA/Concurrency.hpp b/src/libslic3r/SLA/Concurrency.hpp index f82d6a39e..a7b5b6c61 100644 --- a/src/libslic3r/SLA/Concurrency.hpp +++ b/src/libslic3r/SLA/Concurrency.hpp @@ -43,23 +43,36 @@ template<> struct _ccr }); } - template - static T reduce(I from, - I to, - const T & init, - Fn && fn, - MergeFn &&mergefn, - size_t granularity = 1) + template + static T reduce(I from, + I to, + const T &init, + MergeFn &&mergefn, + AccessFn &&access, + size_t granularity = 1 + ) { return tbb::parallel_reduce( tbb::blocked_range{from, to, granularity}, init, [&](const auto &range, T subinit) { T acc = subinit; - loop_(range, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); }); + loop_(range, [&](auto &i) { acc = mergefn(acc, access(i)); }); return acc; }, std::forward(mergefn)); } + + template + static IteratorOnly reduce(I from, + I to, + const T & init, + MergeFn &&mergefn, + size_t granularity = 1) + { + return reduce( + from, to, init, std::forward(mergefn), + [](typename I::value_type &i) { return i; }, granularity); + } }; template<> struct _ccr @@ -92,18 +105,31 @@ public: loop_(from, to, std::forward(fn)); } - template - static IntegerOnly reduce(I from, - I to, - const T & init, - Fn && fn, - MergeFn &&mergefn, - size_t /*granularity*/ = 1) + template + static T reduce(I from, + I to, + const T & init, + MergeFn &&mergefn, + AccessFn &&access, + size_t /*granularity*/ = 1 + ) { T acc = init; - loop_(from, to, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); }); + loop_(from, to, [&](auto &i) { acc = mergefn(acc, access(i)); }); return acc; } + + template + static IteratorOnly reduce(I from, + I to, + const T &init, + MergeFn &&mergefn, + size_t /*granularity*/ = 1 + ) + { + return reduce(from, to, init, std::forward(mergefn), + [](typename I::value_type &i) { return i; }); + } }; using ccr = _ccr; diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index 723e50eeb..b7bed1c91 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -31,7 +31,7 @@ VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) { return vmap; } -// Find transformed mesh ground level without copy and with parallell reduce. +// Find transformed mesh ground level without copy and with parallel reduce. double find_ground_level(const TriangleMesh &mesh, const Transform3d & tr, size_t threads) @@ -40,15 +40,13 @@ double find_ground_level(const TriangleMesh &mesh, auto minfn = [](double a, double b) { return std::min(a, b); }; - auto findminz = [&mesh, &tr] (size_t vi, double submin) { - Vec3d v = tr * mesh.its.vertices[vi].template cast(); - return std::min(submin, v.z()); + auto accessfn = [&mesh, &tr] (size_t vi) { + return (tr * mesh.its.vertices[vi].template cast()).z(); }; double zmin = mesh.its.vertices.front().z(); - - return ccr_par::reduce(size_t(0), vsize, zmin, findminz, minfn, - vsize / threads); + size_t granularity = vsize / threads; + return ccr_par::reduce(size_t(0), vsize, zmin, minfn, accessfn, granularity); } // Try to guess the number of support points needed to support a mesh @@ -65,7 +63,7 @@ double calculate_model_supportedness(const TriangleMesh & mesh, double zmin = find_ground_level(mesh, tr, Nthr); - auto score_mergefn = [&mesh, &tr, zmin](size_t fi, double subscore) { + auto accessfn = [&mesh, &tr, zmin](size_t fi) { static const Vec3d DOWN = {0., 0., -1.}; @@ -83,21 +81,18 @@ double calculate_model_supportedness(const TriangleMesh & mesh, double zlvl = zmin + 0.1; if (p1.z() <= zlvl && p2.z() <= zlvl && p3.z() <= zlvl) { // score += area * POINTS_PER_UNIT_AREA; - return subscore; + return 0.; } double phi = 1. - std::acos(N.dot(DOWN)) / PI; - phi = phi * (phi > 0.5); +// phi = phi * (phi > 0.5); // std::cout << "area: " << area << std::endl; - subscore += area * POINTS_PER_UNIT_AREA * phi; - - return subscore; + return area * POINTS_PER_UNIT_AREA * phi; }; - double score = ccr_seq::reduce(size_t(0), facesize, 0., score_mergefn, - std::plus{}, facesize / Nthr); + double score = ccr_par::reduce(size_t(0), facesize, 0., std::plus{}, accessfn, facesize / Nthr); return score / mesh.its.indices.size(); } @@ -107,7 +102,7 @@ std::array find_best_rotation(const ModelObject& modelobj, std::function statuscb, std::function stopcond) { - static const unsigned MAX_TRIES = 100; + static const unsigned MAX_TRIES = 10000; // return value std::array rot; @@ -158,10 +153,10 @@ std::array find_best_rotation(const ModelObject& modelobj, .max_iterations(max_tries) .rel_score_diff(1e-6) .stop_condition(stopcond), - 10 /*grid size*/); + 100 /*grid size*/); - // We are searching rotations around the three axes x, y, z. Thus the - // problem becomes a 3 dimensional optimization task. + // 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 b = opt::Bound{-PI, PI}; diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt index 30b93eafc..501af0c6f 100644 --- a/tests/libslic3r/CMakeLists.txt +++ b/tests/libslic3r/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable(${_TEST_NAME}_tests test_marchingsquares.cpp test_timeutils.cpp test_voronoi.cpp + test_optimizers.cpp test_png_io.cpp test_timeutils.cpp ) diff --git a/tests/libslic3r/test_optimizers.cpp b/tests/libslic3r/test_optimizers.cpp new file mode 100644 index 000000000..6e84f6a69 --- /dev/null +++ b/tests/libslic3r/test_optimizers.cpp @@ -0,0 +1,59 @@ +#include +#include + +#include + +#include + +void check_opt_result(double score, double ref, double abs_err, double rel_err) +{ + double abs_diff = std::abs(score - ref); + double rel_diff = std::abs(abs_diff / std::abs(ref)); + + bool abs_reached = abs_diff < abs_err; + bool rel_reached = rel_diff < rel_err; + bool precision_reached = abs_reached || rel_reached; + REQUIRE(precision_reached); +} + +template void test_sin(Opt &&opt) +{ + using namespace Slic3r::opt; + + auto optfunc = [](const auto &in) { + auto [phi] = in; + + return std::sin(phi); + }; + + auto init = initvals({PI}); + auto optbounds = bounds({ {0., 2 * PI}}); + + Result result_min = opt.to_min().optimize(optfunc, init, optbounds); + Result result_max = opt.to_max().optimize(optfunc, init, optbounds); + + check_opt_result(result_min.score, -1., 1e-2, 1e-4); + check_opt_result(result_max.score, 1., 1e-2, 1e-4); +} + +template void test_sphere_func(Opt &&opt) +{ + using namespace Slic3r::opt; + + Result result = opt.to_min().optimize([](const auto &in) { + auto [x, y] = in; + + return x * x + y * y + 1.; + }, initvals({.6, -0.2}), bounds({{-1., 1.}, {-1., 1.}})); + + check_opt_result(result.score, 1., 1e-2, 1e-4); +} + +TEST_CASE("Test brute force optimzer for basic 1D and 2D functions", "[Opt]") { + using namespace Slic3r::opt; + + Optimizer opt; + + test_sin(opt); + test_sphere_func(opt); +} diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp index dad2b9097..1575ee0e6 100644 --- a/tests/sla_print/sla_print_tests.cpp +++ b/tests/sla_print/sla_print_tests.cpp @@ -5,6 +5,7 @@ #include "sla_test_utils.hpp" #include +#include namespace { @@ -239,3 +240,14 @@ TEST_CASE("halfcone test", "[halfcone]") { m.require_shared_vertices(); m.WriteOBJFile("Halfcone.obj"); } + +TEST_CASE("Test concurrency") +{ + std::vector vals = grid(0., 100., 10.); + + double ref = std::accumulate(vals.begin(), vals.end(), 0.); + + double s = sla::ccr_par::reduce(vals.begin(), vals.end(), 0., std::plus{}); + + REQUIRE(s == Approx(ref)); +} From b4b9af410037fff27cfa0682e1deeb5744515d86 Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Tue, 1 Sep 2020 20:08:42 +0200 Subject: [PATCH 4/9] cosmethics Comments and cosmethics --- src/libslic3r/CMakeLists.txt | 1 + .../Optimize/BruteforceOptimizer.hpp | 26 ++++++++++++++----- src/libslic3r/SLA/Rotfinder.cpp | 20 +++++--------- 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 263920ecb..e30811133 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -217,6 +217,7 @@ add_library(libslic3r STATIC MarchingSquares.hpp Optimize/Optimizer.hpp Optimize/NLoptOptimizer.hpp + Optimize/BruteforceOptimizer.hpp ${OpenVDBUtils_SOURCES} SLA/Pad.hpp SLA/Pad.cpp diff --git a/src/libslic3r/Optimize/BruteforceOptimizer.hpp b/src/libslic3r/Optimize/BruteforceOptimizer.hpp index 960676e7b..2daef538e 100644 --- a/src/libslic3r/Optimize/BruteforceOptimizer.hpp +++ b/src/libslic3r/Optimize/BruteforceOptimizer.hpp @@ -8,14 +8,18 @@ namespace Slic3r { namespace opt { namespace detail { // Implementing a bruteforce optimizer +// Return the number of iterations needed to reach a specific grid position (idx) template -constexpr long num_iter(const std::array &idx, size_t gridsz) +long num_iter(const std::array &idx, size_t gridsz) { long ret = 0; for (size_t i = 0; i < N; ++i) ret += idx[i] * std::pow(gridsz, i); return ret; } +// Implementation of a grid search where the search interval is sampled in +// equidistant points for each dimension. Grid size determines the number of +// samples for one dimension so the number of function calls is gridsize ^ dimension. struct AlgBurteForce { bool to_min; StopCriteria stc; @@ -23,6 +27,11 @@ struct AlgBurteForce { AlgBurteForce(const StopCriteria &cr, size_t gs): stc{cr}, gridsz{gs} {} + // This function is called recursively for each dimension and generates + // the grid values for the particular dimension. If D is less than zero, + // the object function input values are generated for each dimension and it + // can be evaluated. The current best score is compared with the newly + // returned score and changed appropriately. template bool run(std::array &idx, Result &result, @@ -32,11 +41,12 @@ struct AlgBurteForce { { if (stc.stop_condition()) return false; - if constexpr (D < 0) { + if constexpr (D < 0) { // Let's evaluate fn Input inp; auto max_iter = stc.max_iterations(); - if (max_iter && num_iter(idx, gridsz) >= max_iter) return false; + if (max_iter && num_iter(idx, gridsz) >= max_iter) + return false; for (size_t d = 0; d < N; ++d) { const Bound &b = bounds[d]; @@ -45,12 +55,13 @@ struct AlgBurteForce { } auto score = fn(inp); - if (cmp(score, result.score)) { + if (cmp(score, result.score)) { // Change current score to the new double absdiff = std::abs(score - result.score); result.score = score; result.optimum = inp; + // Check if the required precision is reached. if (absdiff < stc.abs_score_diff() || absdiff < stc.rel_score_diff() * std::abs(score)) return false; @@ -58,9 +69,10 @@ struct AlgBurteForce { } else { for (size_t i = 0; i < gridsz; ++i) { - idx[D] = i; + idx[D] = i; // Mark the current grid position and dig down if (!run(idx, result, bounds, std::forward(fn), - std::forward(cmp))) return false; + std::forward(cmp))) + return false; } } @@ -90,7 +102,7 @@ struct AlgBurteForce { } }; -} // namespace bruteforce_detail +} // namespace detail using AlgBruteForce = detail::AlgBurteForce; diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index b7bed1c91..e8cc7df1b 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -6,14 +6,11 @@ #include #include #include -#include -#include #include "Model.hpp" #include -namespace Slic3r { -namespace sla { +namespace Slic3r { namespace sla { using VertexFaceMap = std::vector>; @@ -51,7 +48,6 @@ double find_ground_level(const TriangleMesh &mesh, // Try to guess the number of support points needed to support a mesh double calculate_model_supportedness(const TriangleMesh & mesh, -// const VertexFaceMap &vmap, const Transform3d & tr) { static constexpr double POINTS_PER_UNIT_AREA = 1.; @@ -111,8 +107,6 @@ std::array find_best_rotation(const ModelObject& modelobj, // rotations TriangleMesh mesh = modelobj.raw_mesh(); mesh.require_shared_vertices(); -// auto vmap = create_vertex_face_map(mesh); -// simplify_mesh(mesh); // For current iteration number unsigned status = 0; @@ -147,21 +141,20 @@ std::array find_best_rotation(const ModelObject& modelobj, return score; }; - // Firing up the genetic optimizer. For now it uses the nlopt library. - + // Firing up the optimizer. opt::Optimizer solver(opt::StopCriteria{} .max_iterations(max_tries) .rel_score_diff(1e-6) .stop_condition(stopcond), - 100 /*grid size*/); + std::sqrt(max_tries)/*grid size*/); // 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 b = opt::Bound{-PI, PI}; - // Now we start the optimization process with initial angles (0, 0, 0) - auto result = solver.to_min().optimize(objfunc, opt::initvals({0.0, 0.0}), + // Now we start the optimization process with initial angles (0, 0) + auto result = solver.to_min().optimize(objfunc, opt::initvals({0., 0.}), opt::bounds({b, b})); // Save the result and fck off @@ -171,5 +164,4 @@ std::array find_best_rotation(const ModelObject& modelobj, return rot; } -} -} +}} // namespace Slic3r::sla From 9f3e7617d86cd5c53a161a606c113b1c2e2b658a Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Tue, 1 Sep 2020 20:08:59 +0200 Subject: [PATCH 5/9] Add Imgui popup for rotation gizmo under SLA --- src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp | 64 +++++++++++++++++++++++++ src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp | 64 ++++++++++++++++++------- src/slic3r/GUI/ImGuiWrapper.cpp | 4 +- 3 files changed, 112 insertions(+), 20 deletions(-) diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp index f3e565686..8365875d9 100644 --- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp +++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp @@ -1,9 +1,15 @@ // Include GLGizmoBase.hpp before I18N.hpp as it includes some libigl code, which overrides our localization "L" macro. #include "GLGizmoRotate.hpp" #include "slic3r/GUI/GLCanvas3D.hpp" +#include "slic3r/GUI/ImGuiWrapper.hpp" #include +#include "slic3r/GUI/GUI_App.hpp" +#include "slic3r/GUI/GUI.hpp" +#include "libslic3r/PresetBundle.hpp" + +#include "libslic3r/SLA/Rotfinder.hpp" namespace Slic3r { namespace GUI { @@ -194,6 +200,64 @@ void GLGizmoRotate::on_render_for_picking() const glsafe(::glPopMatrix()); } +GLGizmoRotate3D::RotoptimzeWindow::RotoptimzeWindow(ImGuiWrapper * imgui, + State & state, + const Alignment &alignment) + : m_imgui{imgui} +{ + imgui->begin(_L("Rotation"), ImGuiWindowFlags_NoMove | + ImGuiWindowFlags_AlwaysAutoResize | + ImGuiWindowFlags_NoCollapse); + + // adjust window position to avoid overlap the view toolbar + float win_h = ImGui::GetWindowHeight(); + float x = alignment.x, y = alignment.y; + y = std::min(y, alignment.bottom_limit - win_h); + ImGui::SetWindowPos(ImVec2(x, y), ImGuiCond_Always); + + ImGui::SliderFloat(_L("Accuracy").c_str(), &state.accuracy, 0.01f, 1.f, "%.1f"); + + if (imgui->button(_L("Optimize orientation"))) { + std::cout << "Blip" << std::endl; + } + + static const std::vector options = { + _L("Least supports").ToStdString(), + _L("Suface quality").ToStdString() + }; + + if (imgui->combo(_L("Choose method"), options, state.method) ) { + std::cout << "method: " << state.method << std::endl; + } + + +} + +GLGizmoRotate3D::RotoptimzeWindow::~RotoptimzeWindow() +{ + m_imgui->end(); +} + +void GLGizmoRotate3D::on_render_input_window(float x, float y, float bottom_limit) +{ + if (wxGetApp().preset_bundle->printers.get_edited_preset().printer_technology() != ptSLA) + return; + +// m_rotoptimizewin_state.mobj = ; + RotoptimzeWindow popup{m_imgui, m_rotoptimizewin_state, {x, y, bottom_limit}}; + +// if ((last_h != win_h) || (last_y != y)) +// { +// // ask canvas for another frame to render the window in the correct position +// m_parent.request_extra_frame(); +// if (last_h != win_h) +// last_h = win_h; +// if (last_y != y) +// last_y = y; +// } + +} + void GLGizmoRotate::render_circle() const { ::glBegin(GL_LINE_LOOP); diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp index 7365a20c3..c547dfbc0 100644 --- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp +++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp @@ -52,12 +52,12 @@ public: std::string get_tooltip() const override; protected: - virtual bool on_init(); - virtual std::string on_get_name() const { return ""; } - virtual void on_start_dragging(); - virtual void on_update(const UpdateData& data); - virtual void on_render() const; - virtual void on_render_for_picking() const; + bool on_init() override; + std::string on_get_name() const override { return ""; } + void on_start_dragging() override; + void on_update(const UpdateData& data) override; + void on_render() const override; + void on_render_for_picking() const override; private: void render_circle() const; @@ -94,46 +94,74 @@ public: } protected: - virtual bool on_init(); - virtual std::string on_get_name() const; - virtual void on_set_state() + bool on_init() override; + std::string on_get_name() const override; + void on_set_state() override { for (GLGizmoRotate& g : m_gizmos) g.set_state(m_state); } - virtual void on_set_hover_id() + void on_set_hover_id() override { for (int i = 0; i < 3; ++i) m_gizmos[i].set_hover_id((m_hover_id == i) ? 0 : -1); } - virtual void on_enable_grabber(unsigned int id) + void on_enable_grabber(unsigned int id) override { if (id < 3) m_gizmos[id].enable_grabber(0); } - virtual void on_disable_grabber(unsigned int id) + void on_disable_grabber(unsigned int id) override { if (id < 3) m_gizmos[id].disable_grabber(0); } - virtual bool on_is_activable() const; - virtual void on_start_dragging(); - virtual void on_stop_dragging(); - virtual void on_update(const UpdateData& data) + bool on_is_activable() const override; + void on_start_dragging() override; + void on_stop_dragging() override; + void on_update(const UpdateData& data) override { for (GLGizmoRotate& g : m_gizmos) { g.update(data); } } - virtual void on_render() const; - virtual void on_render_for_picking() const + void on_render() const override; + void on_render_for_picking() const override { for (const GLGizmoRotate& g : m_gizmos) { g.render_for_picking(); } } + + void on_render_input_window(float x, float y, float bottom_limit) override; +private: + + class RotoptimzeWindow { + ImGuiWrapper *m_imgui = nullptr; + + public: + struct State { + enum Metods { mMinSupportPoints, mLegacy }; + + float accuracy = 1.f; + int method = mMinSupportPoints; + ModelObject *mobj = nullptr; + }; + + struct Alignment { float x, y, bottom_limit; }; + + RotoptimzeWindow(ImGuiWrapper *imgui, State &settings, const Alignment &bottom_limit); + ~RotoptimzeWindow(); + + RotoptimzeWindow(const RotoptimzeWindow&) = delete; + RotoptimzeWindow(RotoptimzeWindow &&) = delete; + RotoptimzeWindow& operator=(const RotoptimzeWindow &) = delete; + RotoptimzeWindow& operator=(RotoptimzeWindow &&) = delete; + }; + + RotoptimzeWindow::State m_rotoptimizewin_state = {}; }; } // namespace GUI diff --git a/src/slic3r/GUI/ImGuiWrapper.cpp b/src/slic3r/GUI/ImGuiWrapper.cpp index e839fdf9b..0fecc822d 100644 --- a/src/slic3r/GUI/ImGuiWrapper.cpp +++ b/src/slic3r/GUI/ImGuiWrapper.cpp @@ -425,10 +425,10 @@ bool ImGuiWrapper::combo(const wxString& label, const std::vector& text(label); ImGui::SameLine(); - int selection_out = -1; + int selection_out = selection; bool res = false; - const char *selection_str = selection < (int)options.size() ? options[selection].c_str() : ""; + const char *selection_str = selection < int(options.size()) && selection >= 0 ? options[selection].c_str() : ""; if (ImGui::BeginCombo("", selection_str)) { for (int i = 0; i < (int)options.size(); i++) { if (ImGui::Selectable(options[i].c_str(), i == selection)) { From 0d4c67b9a34ed5699c2676ebb079c681c763b39a Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Fri, 4 Sep 2020 17:51:22 +0200 Subject: [PATCH 6/9] Mostly working, inefficiencies remain, status indication partly broken --- src/libslic3r/SLA/Rotfinder.cpp | 316 ++++++++++++++++-------- src/libslic3r/SLA/Rotfinder.hpp | 13 +- src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp | 32 +-- src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp | 7 +- src/slic3r/GUI/Jobs/RotoptimizeJob.cpp | 20 +- 5 files changed, 264 insertions(+), 124 deletions(-) diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index e8cc7df1b..e033009aa 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -5,27 +5,23 @@ #include #include #include -#include + +#include "libslic3r/SLAPrint.hpp" +#include "libslic3r/PrintConfig.hpp" + +#include #include "Model.hpp" #include namespace Slic3r { namespace sla { -using VertexFaceMap = std::vector>; +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(); -VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) { - std::vector> vmap(mesh.its.vertices.size()); - - size_t fi = 0; - for (const Vec3i &tri : mesh.its.indices) { - for (int vi = 0; vi < tri.size(); ++vi) { - auto from = vmap[tri(vi)].begin(), to = vmap[tri(vi)].end(); - vmap[tri(vi)].insert(std::lower_bound(from, to, fi), fi); - } - } - - return vmap; + return opt_elevation < EPSILON || opt_padaround; } // Find transformed mesh ground level without copy and with parallel reduce. @@ -41,62 +37,163 @@ double find_ground_level(const TriangleMesh &mesh, return (tr * mesh.its.vertices[vi].template cast()).z(); }; - double zmin = mesh.its.vertices.front().z(); + double zmin = std::numeric_limits::max(); size_t granularity = vsize / threads; return ccr_par::reduce(size_t(0), vsize, zmin, minfn, accessfn, granularity); } -// Try to guess the number of support points needed to support a mesh -double calculate_model_supportedness(const TriangleMesh & mesh, - const Transform3d & tr) +// Get the vertices of a triangle directly in an array of 3 points +std::array get_triangle_vertices(const TriangleMesh &mesh, + size_t faceidx) { - static constexpr double POINTS_PER_UNIT_AREA = 1.; - - if (mesh.its.vertices.empty()) return std::nan(""); - - size_t Nthr = std::thread::hardware_concurrency(); - size_t facesize = mesh.its.indices.size(); - - double zmin = find_ground_level(mesh, tr, Nthr); - - auto accessfn = [&mesh, &tr, zmin](size_t fi) { - - static const Vec3d DOWN = {0., 0., -1.}; - - const auto &face = mesh.its.indices[fi]; - Vec3d p1 = tr * mesh.its.vertices[face(0)].template cast(); - Vec3d p2 = tr * mesh.its.vertices[face(1)].template cast(); - Vec3d p3 = tr * mesh.its.vertices[face(2)].template cast(); - - Vec3d U = p2 - p1; - Vec3d V = p3 - p1; - Vec3d C = U.cross(V); - Vec3d N = C.normalized(); - double area = 0.5 * C.norm(); - - double zlvl = zmin + 0.1; - if (p1.z() <= zlvl && p2.z() <= zlvl && p3.z() <= zlvl) { - // score += area * POINTS_PER_UNIT_AREA; - return 0.; - } - - double phi = 1. - std::acos(N.dot(DOWN)) / PI; -// phi = phi * (phi > 0.5); - - // std::cout << "area: " << area << std::endl; - - return area * POINTS_PER_UNIT_AREA * phi; - }; - - double score = ccr_par::reduce(size_t(0), facesize, 0., std::plus{}, accessfn, facesize / Nthr); - - return score / mesh.its.indices.size(); + const auto &face = mesh.its.indices[faceidx]; + return {Vec3d{mesh.its.vertices[face(0)].cast()}, + Vec3d{mesh.its.vertices[face(1)].cast()}, + Vec3d{mesh.its.vertices[face(2)].cast()}}; } -std::array find_best_rotation(const ModelObject& modelobj, - float accuracy, - std::function statuscb, - std::function stopcond) +std::array get_transformed_triangle(const TriangleMesh &mesh, + const Transform3d & tr, + size_t faceidx) +{ + const auto &tri = get_triangle_vertices(mesh, faceidx); + return {tr * tri[0], tr * tri[1], tr * tri[2]}; +} + +// Get area and normal of a triangle +struct Face { Vec3d normal; double area; }; +inline Face facestats(const std::array &triangle) +{ + Vec3d U = triangle[1] - triangle[0]; + Vec3d V = triangle[2] - triangle[0]; + Vec3d C = U.cross(V); + Vec3d N = C.normalized(); + double area = 0.5 * C.norm(); + + return {N, area}; +} + +inline const Vec3d DOWN = {0., 0., -1.}; +constexpr double POINTS_PER_UNIT_AREA = 1.; + +inline double get_score(const Face &fc) +{ + double phi = 1. - std::acos(fc.normal.dot(DOWN)) / PI; + phi = phi * (phi > 0.5); + phi = phi * phi * phi; + + return fc.area * POINTS_PER_UNIT_AREA * phi; +} + +template +double sum_score(AccessFn &&accessfn, size_t facecount, size_t Nthreads) +{ + double initv = 0.; + auto mergefn = std::plus{}; + size_t grainsize = facecount / Nthreads; + size_t from = 0, to = facecount; + + return ccr_par::reduce(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 Transform3d &tr) +{ + if (mesh.its.vertices.empty()) return std::nan(""); + + auto accessfn = [&mesh, &tr](size_t fi) { + Face fc = facestats(get_transformed_triangle(mesh, tr, fi)); + return get_score(fc); + }; + + size_t facecount = mesh.its.indices.size(); + size_t Nthreads = std::thread::hardware_concurrency(); + return sum_score(accessfn, facecount, Nthreads) / facecount; +} + +double get_model_supportedness_onfloor(const TriangleMesh &mesh, + const Transform3d & tr) +{ + if (mesh.its.vertices.empty()) return std::nan(""); + + size_t Nthreads = std::thread::hardware_concurrency(); + + double zmin = find_ground_level(mesh, tr, Nthreads); + double zlvl = zmin + 0.1; // 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); + Face fc = facestats(tri); + + if (tri[0].z() <= zlvl && tri[1].z() <= zlvl && tri[2].z() <= zlvl) + return -fc.area * POINTS_PER_UNIT_AREA; + + return get_score(fc); + }; + + size_t facecount = mesh.its.indices.size(); + return sum_score(accessfn, facecount, Nthreads) / facecount; +} + +using XYRotation = std::array; + +// prepare the rotation transformation +Transform3d to_transform3d(const XYRotation &rot) +{ + Transform3d rt = Transform3d::Identity(); + rt.rotate(Eigen::AngleAxisd(rot[1], Vec3d::UnitY())); + rt.rotate(Eigen::AngleAxisd(rot[0], Vec3d::UnitX())); + return rt; +} + +XYRotation from_transform3d(const Transform3d &tr) +{ + Vec3d rot3d = Geometry::Transformation {tr}.get_rotation(); + return {rot3d.x(), rot3d.y()}; +} + +// Find the best score from a set of function inputs. Evaluate for every point. +template +std::array find_min_score(Fn &&fn, Cmp &&cmp, It from, It to) +{ + std::array ret; + + double score = std::numeric_limits::max(); + + for (auto it = from; it != to; ++it) { + double sc = fn(*it); + if (cmp(sc, score)) { + score = sc; + ret = *it; + } + } + + return ret; +} + +// collect the rotations for each face of the convex hull +std::vector get_chull_rotations(const TriangleMesh &mesh) +{ + TriangleMesh chull = mesh.convex_hull_3d(); + chull.require_shared_vertices(); + + size_t facecount = chull.its.indices.size(); + auto inputs = reserve_vector(facecount); + + for (size_t fi = 0; fi < facecount; ++fi) { + Face fc = facestats(get_triangle_vertices(chull, fi)); + + auto q = Eigen::Quaterniond{}.FromTwoVectors(fc.normal, DOWN); + inputs.emplace_back(from_transform3d(Transform3d::Identity() * q)); + } + + return inputs; +} + +XYRotation find_best_rotation(const SLAPrintObject & po, + float accuracy, + std::function statuscb, + std::function stopcond) { static const unsigned MAX_TRIES = 10000; @@ -105,10 +202,10 @@ std::array find_best_rotation(const ModelObject& modelobj, // We will use only one instance of this converted mesh to examine different // rotations - TriangleMesh mesh = modelobj.raw_mesh(); + TriangleMesh mesh = po.model_object()->raw_mesh(); mesh.require_shared_vertices(); - // For current iteration number + // To keep track of the number of iterations unsigned status = 0; // The maximum number of iterations @@ -117,51 +214,66 @@ std::array find_best_rotation(const ModelObject& modelobj, // call status callback with zero, because we are at the start statuscb(status); - // So this is the object function which is called by the solver many times - // It has to yield a single value representing the current score. We will - // call the status callback in each iteration but the actual value may be - // the same for subsequent iterations (status goes from 0 to 100 but - // iterations can be many more) - auto objfunc = [&mesh, &status, &statuscb, &stopcond, max_tries] - (const opt::Input<2> &in) - { - std::cout << "in: " << in[0] << " " << in[1] << std::endl; - - // prepare the rotation transformation - Transform3d rt = Transform3d::Identity(); - rt.rotate(Eigen::AngleAxisd(in[1], Vec3d::UnitY())); - rt.rotate(Eigen::AngleAxisd(in[0], Vec3d::UnitX())); - - double score = sla::calculate_model_supportedness(mesh, rt); - std::cout << score << std::endl; - + auto statusfn = [&statuscb, &status, max_tries] { // report status - if(!stopcond()) statuscb( unsigned(++status * 100.0/max_tries) ); - - return score; + statuscb(unsigned(++status * 100.0/max_tries) ); }; - // Firing up the optimizer. - opt::Optimizer solver(opt::StopCriteria{} - .max_iterations(max_tries) - .rel_score_diff(1e-6) - .stop_condition(stopcond), - std::sqrt(max_tries)/*grid size*/); + // Different search methods have to be used depending on the model elevation + if (is_on_floor(po)) { - // 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 b = opt::Bound{-PI, PI}; + // If the model can be placed on the bed directly, we only need to + // check the 3D convex hull face rotations. - // Now we start the optimization process with initial angles (0, 0) - auto result = solver.to_min().optimize(objfunc, opt::initvals({0., 0.}), - opt::bounds({b, b})); + auto inputs = get_chull_rotations(mesh); - // Save the result and fck off - rot[0] = std::get<0>(result.optimum); - rot[1] = std::get<1>(result.optimum); + auto cmpfn = [](double a, double b) { return a < b; }; + auto objfn = [&mesh, &statusfn](const XYRotation &rot) { + statusfn(); + // We actually need the reverserotation to make the object lie on + // this face + Transform3d tr = to_transform3d(rot); + return get_model_supportedness_onfloor(mesh, tr); + }; + + rot = find_min_score<2>(objfn, cmpfn, inputs.begin(), inputs.end()); + } else { + + // Preparing the optimizer. + size_t grid_size = std::sqrt(max_tries); + opt::Optimizer solver(opt::StopCriteria{} + .max_iterations(max_tries) + .stop_condition(stopcond), + grid_size); + + // 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_model_supportedness(mesh, to_transform3d(rot)); + }, opt::initvals({0., 0.}), bounds); + + // Save the result and fck off + rot = result.optimum; + + std::cout << "best score: " << result.score << std::endl; + } return rot; } +double get_model_supportedness(const SLAPrintObject &po, const Transform3d &tr) +{ + TriangleMesh mesh = po.model_object()->raw_mesh(); + mesh.require_shared_vertices(); + + return is_on_floor(po) ? get_model_supportedness_onfloor(mesh, tr) : + get_model_supportedness(mesh, tr); +} + }} // namespace Slic3r::sla diff --git a/src/libslic3r/SLA/Rotfinder.hpp b/src/libslic3r/SLA/Rotfinder.hpp index 583703203..4fa529600 100644 --- a/src/libslic3r/SLA/Rotfinder.hpp +++ b/src/libslic3r/SLA/Rotfinder.hpp @@ -4,9 +4,11 @@ #include #include +#include + namespace Slic3r { -class ModelObject; +class SLAPrintObject; namespace sla { @@ -26,13 +28,16 @@ namespace sla { * @return Returns the rotations around each axis (x, y, z) */ std::array find_best_rotation( - const ModelObject& modelobj, + const SLAPrintObject& modelobj, float accuracy = 1.0f, std::function statuscb = [] (unsigned) {}, std::function stopcond = [] () { return false; } ); -} -} +double get_model_supportedness(const SLAPrintObject &mesh, + const Transform3d & tr); + +} // namespace sla +} // namespace Slic3r #endif // SLAROTFINDER_HPP diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp index 8365875d9..77366c633 100644 --- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp +++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp @@ -200,6 +200,8 @@ void GLGizmoRotate::on_render_for_picking() const glsafe(::glPopMatrix()); } + + GLGizmoRotate3D::RotoptimzeWindow::RotoptimzeWindow(ImGuiWrapper * imgui, State & state, const Alignment &alignment) @@ -215,20 +217,26 @@ GLGizmoRotate3D::RotoptimzeWindow::RotoptimzeWindow(ImGuiWrapper * imgui, y = std::min(y, alignment.bottom_limit - win_h); ImGui::SetWindowPos(ImVec2(x, y), ImGuiCond_Always); - ImGui::SliderFloat(_L("Accuracy").c_str(), &state.accuracy, 0.01f, 1.f, "%.1f"); + static constexpr const char * button_txt = L("Optimize orientation"); + static constexpr const char * slider_txt = L("Accuracy"); - if (imgui->button(_L("Optimize orientation"))) { + float button_width = imgui->calc_text_size(_(button_txt)).x; + ImGui::PushItemWidth(100.); + //if (imgui->button(_(button_txt))) { + if (ImGui::ArrowButton(_(button_txt).c_str(), ImGuiDir_Down)){ std::cout << "Blip" << std::endl; } + ImGui::SliderFloat(_(slider_txt).c_str(), &state.accuracy, 0.01f, 1.f, "%.1f"); + static const std::vector options = { _L("Least supports").ToStdString(), _L("Suface quality").ToStdString() }; - if (imgui->combo(_L("Choose method"), options, state.method) ) { - std::cout << "method: " << state.method << std::endl; - } +// if (imgui->combo(_L("Choose method"), options, state.method) ) { +// std::cout << "method: " << state.method << std::endl; +// } } @@ -243,18 +251,10 @@ void GLGizmoRotate3D::on_render_input_window(float x, float y, float bottom_limi if (wxGetApp().preset_bundle->printers.get_edited_preset().printer_technology() != ptSLA) return; -// m_rotoptimizewin_state.mobj = ; - RotoptimzeWindow popup{m_imgui, m_rotoptimizewin_state, {x, y, bottom_limit}}; +// TODO: -// if ((last_h != win_h) || (last_y != y)) -// { -// // ask canvas for another frame to render the window in the correct position -// m_parent.request_extra_frame(); -// if (last_h != win_h) -// last_h = win_h; -// if (last_y != y) -// last_y = y; -// } +// m_rotoptimizewin_state.mobj = ?; +// RotoptimzeWindow popup{m_imgui, m_rotoptimizewin_state, {x, y, bottom_limit}}; } diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp index c547dfbc0..c418c4b31 100644 --- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp +++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp @@ -136,12 +136,14 @@ protected: } void on_render_input_window(float x, float y, float bottom_limit) override; + private: class RotoptimzeWindow { ImGuiWrapper *m_imgui = nullptr; public: + struct State { enum Metods { mMinSupportPoints, mLegacy }; @@ -152,7 +154,10 @@ private: struct Alignment { float x, y, bottom_limit; }; - RotoptimzeWindow(ImGuiWrapper *imgui, State &settings, const Alignment &bottom_limit); + RotoptimzeWindow(ImGuiWrapper * imgui, + State & state, + const Alignment &bottom_limit); + ~RotoptimzeWindow(); RotoptimzeWindow(const RotoptimzeWindow&) = delete; diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp index 3fd86b13f..91a67cfa2 100644 --- a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp +++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp @@ -4,6 +4,7 @@ #include "libslic3r/SLA/Rotfinder.hpp" #include "libslic3r/MinAreaBoundingBox.hpp" #include "libslic3r/Model.hpp" +#include "libslic3r/SLAPrint.hpp" #include "slic3r/GUI/Plater.hpp" @@ -15,9 +16,26 @@ void RotoptimizeJob::process() if (obj_idx < 0) { return; } ModelObject *o = m_plater->model().objects[size_t(obj_idx)]; + const SLAPrintObject *po = m_plater->sla_print().objects()[size_t(obj_idx)]; + + if (!o || !po) return; + + TriangleMesh mesh = o->raw_mesh(); + mesh.require_shared_vertices(); + +// for (auto inst : o->instances) { +// Transform3d tr = Transform3d::Identity(); +// tr.rotate(Eigen::AngleAxisd(inst->get_rotation(Z), Vec3d::UnitZ())); +// tr.rotate(Eigen::AngleAxisd(inst->get_rotation(Y), Vec3d::UnitY())); +// tr.rotate(Eigen::AngleAxisd(inst->get_rotation(X), Vec3d::UnitX())); + +// double score = sla::get_model_supportedness(*po, tr); + +// std::cout << "Model supportedness before: " << score << std::endl; +// } auto r = sla::find_best_rotation( - *o, + *po, 1.f, [this](unsigned s) { if (s < 100) From 3b7ea5587e5b6c7d6938c746dd2855b95d242aa7 Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Fri, 4 Sep 2020 19:43:07 +0200 Subject: [PATCH 7/9] Fix build on win --- src/libslic3r/SLA/Concurrency.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/libslic3r/SLA/Concurrency.hpp b/src/libslic3r/SLA/Concurrency.hpp index a7b5b6c61..300024c76 100644 --- a/src/libslic3r/SLA/Concurrency.hpp +++ b/src/libslic3r/SLA/Concurrency.hpp @@ -5,7 +5,10 @@ #include #include #include + #include +#include + #include namespace Slic3r { From d5271220464fb8da07bdc805c68318ece201d3bb Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Fri, 4 Sep 2020 20:20:06 +0200 Subject: [PATCH 8/9] Performance optimizations and bugfix --- src/libslic3r/SLA/Rotfinder.cpp | 16 ++++++++++++++-- src/slic3r/GUI/Jobs/RotoptimizeJob.cpp | 3 ++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index e033009aa..db8c0b9a8 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -76,12 +76,20 @@ inline Face facestats(const std::array &triangle) inline const Vec3d DOWN = {0., 0., -1.}; constexpr double POINTS_PER_UNIT_AREA = 1.; +// The score function for a particular face inline double get_score(const Face &fc) { + // Simply get the angle (acos of dot product) between the face normal and + // the DOWN vector. double phi = 1. - std::acos(fc.normal.dot(DOWN)) / 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; } @@ -176,6 +184,8 @@ std::vector get_chull_rotations(const TriangleMesh &mesh) { 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(); auto inputs = reserve_vector(facecount); @@ -183,8 +193,10 @@ std::vector get_chull_rotations(const TriangleMesh &mesh) for (size_t fi = 0; fi < facecount; ++fi) { Face fc = facestats(get_triangle_vertices(chull, fi)); - auto q = Eigen::Quaterniond{}.FromTwoVectors(fc.normal, DOWN); - inputs.emplace_back(from_transform3d(Transform3d::Identity() * q)); + if (fc.area > area_threshold) { + auto q = Eigen::Quaterniond{}.FromTwoVectors(fc.normal, DOWN); + inputs.emplace_back(from_transform3d(Transform3d::Identity() * q)); + } } return inputs; diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp index 91a67cfa2..10c09275c 100644 --- a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp +++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp @@ -13,7 +13,8 @@ namespace Slic3r { namespace GUI { void RotoptimizeJob::process() { int obj_idx = m_plater->get_selected_object_idx(); - if (obj_idx < 0) { return; } + if (obj_idx < 0 || m_plater->sla_print().objects().size() <= obj_idx) + return; ModelObject *o = m_plater->model().objects[size_t(obj_idx)]; const SLAPrintObject *po = m_plater->sla_print().objects()[size_t(obj_idx)]; From 20bd7b99f9c8d4ded94e6c01450ffdbfadc36c0a Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Thu, 10 Sep 2020 19:35:26 +0200 Subject: [PATCH 9/9] Significant performance improvements for elevated and non-elevated case Apply bruteforce for elevated models --- src/libslic3r/Fill/FillAdaptive.hpp | 1 + src/libslic3r/SLA/Rotfinder.cpp | 135 +++++++++++++++---------- src/libslic3r/SLA/Rotfinder.hpp | 2 +- src/slic3r/GUI/Jobs/RotoptimizeJob.cpp | 11 +- 4 files changed, 89 insertions(+), 60 deletions(-) diff --git a/src/libslic3r/Fill/FillAdaptive.hpp b/src/libslic3r/Fill/FillAdaptive.hpp index dd7394384..23786530e 100644 --- a/src/libslic3r/Fill/FillAdaptive.hpp +++ b/src/libslic3r/Fill/FillAdaptive.hpp @@ -4,6 +4,7 @@ #include "../AABBTreeIndirect.hpp" #include "FillBase.hpp" +#include "TriangleMesh.hpp" namespace Slic3r { diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp index db8c0b9a8..937897766 100644 --- a/src/libslic3r/SLA/Rotfinder.cpp +++ b/src/libslic3r/SLA/Rotfinder.cpp @@ -1,11 +1,10 @@ #include -#include -//#include -#include #include #include +#include + #include "libslic3r/SLAPrint.hpp" #include "libslic3r/PrintConfig.hpp" @@ -61,23 +60,25 @@ std::array get_transformed_triangle(const TriangleMesh &mesh, } // Get area and normal of a triangle -struct Face { Vec3d normal; double area; }; -inline Face facestats(const std::array &triangle) -{ - Vec3d U = triangle[1] - triangle[0]; - Vec3d V = triangle[2] - triangle[0]; - Vec3d C = U.cross(V); - Vec3d N = C.normalized(); - double area = 0.5 * C.norm(); +struct Facestats { + Vec3d normal; + double area; - return {N, area}; -} + explicit Facestats(const std::array &triangle) + { + Vec3d U = triangle[1] - triangle[0]; + Vec3d V = triangle[2] - triangle[0]; + Vec3d C = U.cross(V); + normal = C.normalized(); + area = 0.5 * C.norm(); + } +}; inline const Vec3d DOWN = {0., 0., -1.}; constexpr double POINTS_PER_UNIT_AREA = 1.; // The score function for a particular face -inline double get_score(const Face &fc) +inline double get_score(const Facestats &fc) { // Simply get the angle (acos of dot product) between the face normal and // the DOWN vector. @@ -110,7 +111,7 @@ double get_model_supportedness(const TriangleMesh &mesh, const Transform3d &tr) if (mesh.its.vertices.empty()) return std::nan(""); auto accessfn = [&mesh, &tr](size_t fi) { - Face fc = facestats(get_transformed_triangle(mesh, tr, fi)); + Facestats fc{get_transformed_triangle(mesh, tr, fi)}; return get_score(fc); }; @@ -131,7 +132,7 @@ double get_model_supportedness_onfloor(const TriangleMesh &mesh, auto accessfn = [&mesh, &tr, zlvl](size_t fi) { std::array tri = get_transformed_triangle(mesh, tr, fi); - Face fc = facestats(tri); + Facestats fc{tri}; if (tri[0].z() <= zlvl && tri[1].z() <= zlvl && tri[2].z() <= zlvl) return -fc.area * POINTS_PER_UNIT_AREA; @@ -161,56 +162,91 @@ XYRotation from_transform3d(const Transform3d &tr) } // Find the best score from a set of function inputs. Evaluate for every point. -template -std::array find_min_score(Fn &&fn, Cmp &&cmp, It from, It to) +template +std::array find_min_score(Fn &&fn, It from, It to, StopCond &&stopfn) { std::array ret; double score = std::numeric_limits::max(); - for (auto it = from; it != to; ++it) { - double sc = fn(*it); - if (cmp(sc, score)) { - score = sc; - ret = *it; - } - } + size_t Nthreads = std::thread::hardware_concurrency(); + size_t dist = std::distance(from, to); + std::vector scores(dist, score); + + ccr_par::for_each(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; } // collect the rotations for each face of the convex hull -std::vector get_chull_rotations(const TriangleMesh &mesh) +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.)); + double area_threshold = chull2d_area / (scaled(1e3) * scaled(1.)); size_t facecount = chull.its.indices.size(); - auto inputs = reserve_vector(facecount); + + 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) { - Face fc = facestats(get_triangle_vertices(chull, fi)); + Facestats fc{get_triangle_vertices(chull, fi)}; if (fc.area > area_threshold) { auto q = Eigen::Quaterniond{}.FromTwoVectors(fc.normal, DOWN); - inputs.emplace_back(from_transform3d(Transform3d::Identity() * q)); + XYRotation rot = from_transform3d(Transform3d::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); } } - return inputs; + 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; } -XYRotation find_best_rotation(const SLAPrintObject & po, - float accuracy, - std::function statuscb, - std::function stopcond) +Vec2d find_best_rotation(const SLAPrintObject & po, + float accuracy, + std::function statuscb, + std::function stopcond) { - static const unsigned MAX_TRIES = 10000; + static const unsigned MAX_TRIES = 1000; // return value - std::array rot; + XYRotation rot; // We will use only one instance of this converted mesh to examine different // rotations @@ -226,7 +262,7 @@ XYRotation find_best_rotation(const SLAPrintObject & po, // call status callback with zero, because we are at the start statuscb(status); - auto statusfn = [&statuscb, &status, max_tries] { + auto statusfn = [&statuscb, &status, &max_tries] { // report status statuscb(unsigned(++status * 100.0/max_tries) ); }; @@ -234,29 +270,26 @@ XYRotation find_best_rotation(const SLAPrintObject & po, // 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 inputs = get_chull_rotations(mesh); - - auto cmpfn = [](double a, double b) { return a < b; }; auto objfn = [&mesh, &statusfn](const XYRotation &rot) { statusfn(); - // We actually need the reverserotation to make the object lie on - // this face Transform3d tr = to_transform3d(rot); return get_model_supportedness_onfloor(mesh, tr); }; - rot = find_min_score<2>(objfn, cmpfn, inputs.begin(), inputs.end()); + rot = find_min_score<2>(objfn, inputs.begin(), inputs.end(), stopcond); } else { - // Preparing the optimizer. - size_t grid_size = std::sqrt(max_tries); + 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), - grid_size); + .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. @@ -272,11 +305,9 @@ XYRotation find_best_rotation(const SLAPrintObject & po, // Save the result and fck off rot = result.optimum; - - std::cout << "best score: " << result.score << std::endl; } - return rot; + return {rot[0], rot[1]}; } double get_model_supportedness(const SLAPrintObject &po, const Transform3d &tr) diff --git a/src/libslic3r/SLA/Rotfinder.hpp b/src/libslic3r/SLA/Rotfinder.hpp index 4fa529600..96561a890 100644 --- a/src/libslic3r/SLA/Rotfinder.hpp +++ b/src/libslic3r/SLA/Rotfinder.hpp @@ -27,7 +27,7 @@ namespace sla { * * @return Returns the rotations around each axis (x, y, z) */ -std::array find_best_rotation( +Vec2d find_best_rotation( const SLAPrintObject& modelobj, float accuracy = 1.0f, std::function statuscb = [] (unsigned) {}, diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp index 10c09275c..978ccf8fc 100644 --- a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp +++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp @@ -13,7 +13,7 @@ namespace Slic3r { namespace GUI { void RotoptimizeJob::process() { int obj_idx = m_plater->get_selected_object_idx(); - if (obj_idx < 0 || m_plater->sla_print().objects().size() <= obj_idx) + if (obj_idx < 0 || int(m_plater->sla_print().objects().size()) <= obj_idx) return; ModelObject *o = m_plater->model().objects[size_t(obj_idx)]; @@ -35,15 +35,12 @@ void RotoptimizeJob::process() // std::cout << "Model supportedness before: " << score << std::endl; // } - auto r = sla::find_best_rotation( - *po, - 1.f, + Vec2d r = sla::find_best_rotation(*po, 0.75f, [this](unsigned s) { if (s < 100) - update_status(int(s), - _(L("Searching for optimal orientation"))); + update_status(int(s), _(L("Searching for optimal orientation"))); }, - [this]() { return was_canceled(); }); + [this] () { return was_canceled(); }); double mindist = 6.0; // FIXME