From b4e30cc8ad08de86fee89b947c35cd401ca9021a Mon Sep 17 00:00:00 2001 From: tamasmeszaros Date: Wed, 26 Aug 2020 10:25:09 +0200 Subject: [PATCH] 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);