From 927b81ea9710266e2effb050ca27cb8762239870 Mon Sep 17 00:00:00 2001 From: tamasmeszaros <meszaros.q@gmail.com> Date: Thu, 16 Jul 2020 11:49:30 +0200 Subject: [PATCH] Working small-to-normal support merging Fixed fatal bug with anchors for mini supports Make the optimization cleaner in support generatior Much better widening behaviour Add an optimizer interface and the NLopt implementation into libslic3r New optimizer based only on nlopt C interfase Fix build and tests --- src/libslic3r/CMakeLists.txt | 1 + src/libslic3r/Optimizer.hpp | 369 ++++++++++ src/libslic3r/SLA/IndexedMesh.hpp | 2 +- src/libslic3r/SLA/SupportTreeBuilder.cpp | 5 + src/libslic3r/SLA/SupportTreeBuilder.hpp | 45 +- src/libslic3r/SLA/SupportTreeBuildsteps.cpp | 758 ++++++++++---------- src/libslic3r/SLA/SupportTreeBuildsteps.hpp | 17 +- src/libslic3r/SLA/SupportTreeMesher.hpp | 18 + tests/sla_print/sla_print_tests.cpp | 11 + 9 files changed, 831 insertions(+), 395 deletions(-) create mode 100644 src/libslic3r/Optimizer.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 91da5df5d..58b74402e 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -203,6 +203,7 @@ add_library(libslic3r STATIC SimplifyMeshImpl.hpp SimplifyMesh.cpp MarchingSquares.hpp + Optimizer.hpp ${OpenVDBUtils_SOURCES} SLA/Pad.hpp SLA/Pad.cpp diff --git a/src/libslic3r/Optimizer.hpp b/src/libslic3r/Optimizer.hpp new file mode 100644 index 000000000..dc70abe33 --- /dev/null +++ b/src/libslic3r/Optimizer.hpp @@ -0,0 +1,369 @@ +#ifndef NLOPTOPTIMIZER_HPP +#define NLOPTOPTIMIZER_HPP + +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable: 4244) +#pragma warning(disable: 4267) +#endif +#include <nlopt.h> +#ifdef _MSC_VER +#pragma warning(pop) +#endif + +#include <utility> +#include <tuple> +#include <array> +#include <cmath> +#include <functional> +#include <limits> +#include <cassert> + +namespace Slic3r { namespace opt { + +// A type to hold the complete result of the optimization. +template<size_t N> struct Result { + int resultcode; + std::array<double, N> 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<double>::min(), + double max = std::numeric_limits<double>::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<size_t N> using Input = std::array<double, N>; +template<size_t N> using Bounds = std::array<Bound, N>; + +// 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<bool()> 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<class Fn> StopCriteria & stop_condition(Fn &&cond) + { + m_stop_condition = cond; return *this; + } + + bool stop_condition() { return m_stop_condition(); } +}; + +// Helper to be used in static_assert. +template<class T> struct always_false { enum { value = false }; }; + +// Basic interface to optimizer object +template<class Method, class Enable = void> class Optimizer { +public: + + Optimizer(const StopCriteria &) + { + static_assert(always_false<Method>::value, + "Optimizer unimplemented for given method!"); + } + + Optimizer<Method, Enable> &to_min() { return *this; } + Optimizer<Method, Enable> &to_max() { return *this; } + Optimizer<Method, Enable> &set_criteria(const StopCriteria &) { return *this; } + StopCriteria get_criteria() const { return {}; }; + + template<class Func, size_t N> + Result<N> optimize(Func&& func, + const Input<N> &initvals, + const Bounds<N>& bounds) { return {}; } + + // optional for randomized methods: + void seed(long /*s*/) {} +}; + +namespace detail { + +// Helper types for NLopt algorithm selection in template contexts +template<nlopt_algorithm alg> struct NLoptAlg {}; + +// NLopt can combine multiple algorithms if one is global an other is a local +// method. This is how template specializations can be informed about this fact. +template<nlopt_algorithm gl_alg, nlopt_algorithm lc_alg = NLOPT_LN_NELDERMEAD> +struct NLoptAlgComb {}; + +template<class M> struct IsNLoptAlg { + static const constexpr bool value = false; +}; + +template<nlopt_algorithm a> struct IsNLoptAlg<NLoptAlg<a>> { + static const constexpr bool value = true; +}; + +template<nlopt_algorithm a1, nlopt_algorithm a2> +struct IsNLoptAlg<NLoptAlgComb<a1, a2>> { + static const constexpr bool value = true; +}; + +template<class M, class T = void> +using NLoptOnly = std::enable_if_t<IsNLoptAlg<M>::value, T>; + +// Convert any collection to tuple. This is useful for object functions taking +// an argument list of doubles. Make things cleaner on the call site of +// optimize(). +template<size_t I, std::size_t N, class T, class C> struct to_tuple_ { + static auto call(const C &c) + { + return std::tuple_cat(std::tuple<T>(c[N-I]), + to_tuple_<I-1, N, T, C>::call(c)); + } +}; + +template<size_t N, class T, class C> struct to_tuple_<0, N, T, C> { + static auto call(const C &c) { return std::tuple<>(); } +}; + +// C array to tuple +template<std::size_t N, class T> auto carray_tuple(const T *v) +{ + return to_tuple_<N, N, T, const T*>::call(v); +} + +// Helper to convert C style array to std::array +template<size_t N, class T> auto to_arr(const T (&a) [N]) +{ + std::array<T, N> r; + std::copy(std::begin(a), std::end(a), std::begin(r)); + return r; +} + +enum class OptDir { MIN, MAX }; // Where to optimize + +struct NLopt { // Helper RAII class for nlopt_opt + nlopt_opt ptr = nullptr; + + template<class...A> explicit NLopt(A&&...a) + { + ptr = nlopt_create(std::forward<A>(a)...); + } + + NLopt(const NLopt&) = delete; + NLopt(NLopt&&) = delete; + NLopt& operator=(const NLopt&) = delete; + NLopt& operator=(NLopt&&) = delete; + + ~NLopt() { nlopt_destroy(ptr); } +}; + +template<class Method> class NLoptOpt {}; + +// Optimizers based on NLopt. +template<nlopt_algorithm alg> class NLoptOpt<NLoptAlg<alg>> { +protected: + StopCriteria m_stopcr; + OptDir m_dir; + + template<class Fn> using TOptData = + std::tuple<std::remove_reference_t<Fn>*, NLoptOpt*, nlopt_opt>; + + template<class Fn, size_t N> + static double optfunc(unsigned n, const double *params, + double *gradient, + void *data) + { + assert(n >= N); + + auto tdata = static_cast<TOptData<Fn>*>(data); + + if (std::get<1>(*tdata)->m_stopcr.stop_condition()) + nlopt_force_stop(std::get<2>(*tdata)); + + auto fnptr = std::get<0>(*tdata); + auto funval = carray_tuple<N>(params); + + return std::apply(*fnptr, funval); + } + + template<size_t N> + void set_up(NLopt &nl, const Bounds<N>& bounds) + { + std::array<double, N> lb, ub; + + for (size_t i = 0; i < N; ++i) { + lb[i] = bounds[i].min(); + ub[i] = bounds[i].max(); + } + + nlopt_set_lower_bounds(nl.ptr, lb.data()); + nlopt_set_upper_bounds(nl.ptr, ub.data()); + + double abs_diff = m_stopcr.abs_score_diff(); + double rel_diff = m_stopcr.rel_score_diff(); + double stopval = m_stopcr.stop_score(); + if(!std::isnan(abs_diff)) nlopt_set_ftol_abs(nl.ptr, abs_diff); + if(!std::isnan(rel_diff)) nlopt_set_ftol_rel(nl.ptr, rel_diff); + if(!std::isnan(stopval)) nlopt_set_stopval(nl.ptr, stopval); + + if(this->m_stopcr.max_iterations() > 0) + nlopt_set_maxeval(nl.ptr, this->m_stopcr.max_iterations()); + } + + template<class Fn, size_t N> + Result<N> optimize(NLopt &nl, Fn &&fn, const Input<N> &initvals) + { + Result<N> r; + + TOptData<Fn> data = std::make_tuple(&fn, this, nl.ptr); + + switch(m_dir) { + case OptDir::MIN: + nlopt_set_min_objective(nl.ptr, optfunc<Fn, N>, &data); break; + case OptDir::MAX: + nlopt_set_max_objective(nl.ptr, optfunc<Fn, N>, &data); break; + } + + r.optimum = initvals; + r.resultcode = nlopt_optimize(nl.ptr, r.optimum.data(), &r.score); + + return r; + } + +public: + + template<class Func, size_t N> + Result<N> optimize(Func&& func, + const Input<N> &initvals, + const Bounds<N>& bounds) + { + NLopt nl{alg, N}; + set_up(nl, bounds); + + return optimize(nl, std::forward<Func>(func), initvals); + } + + explicit NLoptOpt(StopCriteria stopcr = {}) : m_stopcr(stopcr) {} + + void set_criteria(const StopCriteria &cr) { m_stopcr = cr; } + const StopCriteria &get_criteria() const noexcept { return m_stopcr; } + void set_dir(OptDir dir) noexcept { m_dir = dir; } + + void seed(long s) { nlopt_srand(s); } +}; + +template<nlopt_algorithm glob, nlopt_algorithm loc> +class NLoptOpt<NLoptAlgComb<glob, loc>>: public NLoptOpt<NLoptAlg<glob>> +{ + using Base = NLoptOpt<NLoptAlg<glob>>; +public: + + template<class Fn, size_t N> + Result<N> optimize(Fn&& f, + const Input<N> &initvals, + const Bounds<N>& bounds) + { + NLopt nl_glob{glob, N}, nl_loc{loc, N}; + + Base::set_up(nl_glob, bounds); + Base::set_up(nl_loc, bounds); + nlopt_set_local_optimizer(nl_glob.ptr, nl_loc.ptr); + + return Base::optimize(nl_glob, std::forward<Fn>(f), initvals); + } + + explicit NLoptOpt(StopCriteria stopcr = {}) : Base{stopcr} {} +}; + +} // namespace detail; + +// Optimizers based on NLopt. +template<class M> class Optimizer<M, detail::NLoptOnly<M>> { + detail::NLoptOpt<M> m_opt; + +public: + + Optimizer& to_max() { m_opt.set_dir(detail::OptDir::MAX); return *this; } + Optimizer& to_min() { m_opt.set_dir(detail::OptDir::MIN); return *this; } + + template<class Func, size_t N> + Result<N> optimize(Func&& func, + const Input<N> &initvals, + const Bounds<N>& bounds) + { + return m_opt.optimize(std::forward<Func>(func), initvals, bounds); + } + + explicit Optimizer(StopCriteria stopcr = {}) : m_opt(stopcr) {} + + Optimizer &set_criteria(const StopCriteria &cr) + { + m_opt.set_criteria(cr); return *this; + } + + const StopCriteria &get_criteria() const { return m_opt.get_criteria(); } + + void seed(long s) { m_opt.seed(s); } +}; + +template<size_t N> Bounds<N> bounds(const Bound (&b) [N]) { return detail::to_arr(b); } +template<size_t N> Input<N> initvals(const double (&a) [N]) { return detail::to_arr(a); } + +// Predefinded NLopt algorithms that are used in the codebase +using AlgNLoptGenetic = detail::NLoptAlgComb<NLOPT_GN_ESCH>; +using AlgNLoptSubplex = detail::NLoptAlg<NLOPT_LN_SBPLX>; +using AlgNLoptSimplex = detail::NLoptAlg<NLOPT_LN_NELDERMEAD>; + +// Helper defs for pre-crafted global and local optimizers that work well. +using DefaultGlobalOptimizer = Optimizer<AlgNLoptGenetic>; +using DefaultLocalOptimizer = Optimizer<AlgNLoptSubplex>; + +}} // namespace Slic3r::opt + +#endif // NLOPTOPTIMIZER_HPP diff --git a/src/libslic3r/SLA/IndexedMesh.hpp b/src/libslic3r/SLA/IndexedMesh.hpp index b0970608e..a72492b34 100644 --- a/src/libslic3r/SLA/IndexedMesh.hpp +++ b/src/libslic3r/SLA/IndexedMesh.hpp @@ -87,7 +87,7 @@ public: inline Vec3d position() const { return m_source + m_dir * m_t; } inline int face() const { return m_face_id; } inline bool is_valid() const { return m_mesh != nullptr; } - inline bool is_hit() const { return !std::isinf(m_t); } + inline bool is_hit() const { return m_face_id >= 0 && !std::isinf(m_t); } inline const Vec3d& normal() const { assert(is_valid()); diff --git a/src/libslic3r/SLA/SupportTreeBuilder.cpp b/src/libslic3r/SLA/SupportTreeBuilder.cpp index 959093623..daa01ef24 100644 --- a/src/libslic3r/SLA/SupportTreeBuilder.cpp +++ b/src/libslic3r/SLA/SupportTreeBuilder.cpp @@ -156,6 +156,11 @@ const TriangleMesh &SupportTreeBuilder::merged_mesh(size_t steps) const merged.merge(get_mesh(bs, steps)); } + for (auto &bs : m_diffbridges) { + if (ctl().stopcondition()) break; + merged.merge(get_mesh(bs, steps)); + } + for (auto &anch : m_anchors) { if (ctl().stopcondition()) break; merged.merge(get_mesh(anch, steps)); diff --git a/src/libslic3r/SLA/SupportTreeBuilder.hpp b/src/libslic3r/SLA/SupportTreeBuilder.hpp index aa8a4ea83..f29263ca3 100644 --- a/src/libslic3r/SLA/SupportTreeBuilder.hpp +++ b/src/libslic3r/SLA/SupportTreeBuilder.hpp @@ -177,6 +177,14 @@ struct Bridge: public SupportTreeNode { Vec3d get_dir() const { return (endp - startp).normalized(); } }; +struct DiffBridge: public Bridge { + double end_r; + + DiffBridge(const Vec3d &p_s, const Vec3d &p_e, double r_s, double r_e) + : Bridge{p_s, p_e, r_s}, end_r{r_e} + {} +}; + // A wrapper struct around the pad struct Pad { TriangleMesh tmesh; @@ -210,14 +218,15 @@ struct Pad { // merged mesh. It can be retrieved using a dedicated method (pad()) class SupportTreeBuilder: public SupportTree { // For heads it is beneficial to use the same IDs as for the support points. - std::vector<Head> m_heads; - std::vector<size_t> m_head_indices; - std::vector<Pillar> m_pillars; - std::vector<Junction> m_junctions; - std::vector<Bridge> m_bridges; - std::vector<Bridge> m_crossbridges; - std::vector<Pedestal> m_pedestals; - std::vector<Anchor> m_anchors; + std::vector<Head> m_heads; + std::vector<size_t> m_head_indices; + std::vector<Pillar> m_pillars; + std::vector<Junction> m_junctions; + std::vector<Bridge> m_bridges; + std::vector<Bridge> m_crossbridges; + std::vector<DiffBridge> m_diffbridges; + std::vector<Pedestal> m_pedestals; + std::vector<Anchor> m_anchors; Pad m_pad; @@ -228,8 +237,8 @@ class SupportTreeBuilder: public SupportTree { mutable bool m_meshcache_valid = false; mutable double m_model_height = 0; // the full height of the model - template<class...Args> - const Bridge& _add_bridge(std::vector<Bridge> &br, Args&&... args) + template<class BridgeT, class...Args> + const BridgeT& _add_bridge(std::vector<BridgeT> &br, Args&&... args) { std::lock_guard<Mutex> lk(m_mutex); br.emplace_back(std::forward<Args>(args)...); @@ -331,17 +340,6 @@ public: return pillar.id; } - const Pillar& head_pillar(unsigned headid) const - { - std::lock_guard<Mutex> lk(m_mutex); - assert(headid < m_head_indices.size()); - - const Head& h = m_heads[m_head_indices[headid]]; - assert(h.pillar_id >= 0 && h.pillar_id < long(m_pillars.size())); - - return m_pillars[size_t(h.pillar_id)]; - } - template<class...Args> const Junction& add_junction(Args&&... args) { std::lock_guard<Mutex> lk(m_mutex); @@ -374,6 +372,11 @@ public: { return _add_bridge(m_crossbridges, std::forward<Args>(args)...); } + + template<class...Args> const DiffBridge& add_diffbridge(Args&&... args) + { + return _add_bridge(m_diffbridges, std::forward<Args>(args)...); + } Head &head(unsigned id) { diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp index 7f6c034dd..7ed410802 100644 --- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp +++ b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp @@ -1,18 +1,25 @@ #include <libslic3r/SLA/SupportTreeBuildsteps.hpp> #include <libslic3r/SLA/SpatIndex.hpp> -#include <libnest2d/optimizers/nlopt/genetic.hpp> -#include <libnest2d/optimizers/nlopt/subplex.hpp> +#include <libslic3r/Optimizer.hpp> #include <boost/log/trivial.hpp> namespace Slic3r { namespace sla { -using libnest2d::opt::initvals; -using libnest2d::opt::bound; -using libnest2d::opt::StopCriteria; -using libnest2d::opt::GeneticOptimizer; -using libnest2d::opt::SubplexOptimizer; +using Slic3r::opt::initvals; +using Slic3r::opt::bounds; +using Slic3r::opt::StopCriteria; +using Slic3r::opt::Optimizer; +using Slic3r::opt::AlgNLoptSubplex; +using Slic3r::opt::AlgNLoptGenetic; + +StopCriteria get_criteria(const SupportTreeConfig &cfg) +{ + return StopCriteria{} + .rel_score_diff(cfg.optimizer_rel_score_diff) + .max_iterations(cfg.optimizer_max_iterations); +} template<class C, class Hit = IndexedMesh::hit_result> static Hit min_hit(const C &hits) @@ -37,7 +44,7 @@ SupportTreeBuildsteps::SupportTreeBuildsteps(SupportTreeBuilder & builder, { // Prepare the support points in Eigen/IGL format as well, we will use // it mostly in this form. - + long i = 0; for (const SupportPoint &sp : m_support_pts) { m_points.row(i)(X) = double(sp.pos(X)); @@ -51,7 +58,7 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, const SupportableMesh &sm) { if(sm.pts.empty()) return false; - + builder.ground_level = sm.emesh.ground_level() - sm.cfg.object_elevation_mm; SupportTreeBuildsteps alg(builder, sm); @@ -72,46 +79,46 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, NUM_STEPS //... }; - + // Collect the algorithm steps into a nice sequence std::array<std::function<void()>, NUM_STEPS> program = { [] () { // Begin... // Potentially clear up the shared data (not needed for now) }, - + std::bind(&SupportTreeBuildsteps::filter, &alg), - + std::bind(&SupportTreeBuildsteps::add_pinheads, &alg), - + std::bind(&SupportTreeBuildsteps::classify, &alg), - + std::bind(&SupportTreeBuildsteps::routing_to_ground, &alg), - + std::bind(&SupportTreeBuildsteps::routing_to_model, &alg), - + std::bind(&SupportTreeBuildsteps::interconnect_pillars, &alg), - + std::bind(&SupportTreeBuildsteps::merge_result, &alg), - + [] () { // Done }, - + [] () { // Abort } }; - + Steps pc = BEGIN; - + if(sm.cfg.ground_facing_only) { program[ROUTING_NONGROUND] = []() { BOOST_LOG_TRIVIAL(info) << "Skipping model-facing supports as requested."; }; } - + // Let's define a simple automaton that will run our program. auto progress = [&builder, &pc] () { static const std::array<std::string, NUM_STEPS> stepstr { @@ -126,7 +133,7 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, "Done", "Abort" }; - + static const std::array<unsigned, NUM_STEPS> stepstate { 0, 10, @@ -139,9 +146,9 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, 100, 0 }; - + if(builder.ctl().stopcondition()) pc = ABORT; - + switch(pc) { case BEGIN: pc = FILTER; break; case FILTER: pc = PINHEADS; break; @@ -155,16 +162,16 @@ bool SupportTreeBuildsteps::execute(SupportTreeBuilder & builder, case ABORT: break; default: ; } - + builder.ctl().statuscb(stepstate[pc], stepstr[pc]); }; - + // Just here we run the computation... while(pc < DONE) { progress(); program[pc](); } - + return pc == ABORT; } @@ -177,48 +184,48 @@ IndexedMesh::hit_result SupportTreeBuildsteps::pinhead_mesh_intersect( double sd) { static const size_t SAMPLES = 8; - + // Move away slightly from the touching point to avoid raycasting on the // inner surface of the mesh. - + auto& m = m_mesh; using HitResult = IndexedMesh::hit_result; - + // Hit results std::array<HitResult, SAMPLES> hits; - + struct Rings { double rpin; double rback; Vec3d spin; Vec3d sback; PointRing<SAMPLES> ring; - + Vec3d backring(size_t idx) { return ring.get(idx, sback, rback); } Vec3d pinring(size_t idx) { return ring.get(idx, spin, rpin); } } rings {r_pin + sd, r_back + sd, s, s + width * dir, dir}; - + // We will shoot multiple rays from the head pinpoint in the direction // of the pinhead robe (side) surface. The result will be the smallest // hit distance. - - ccr::enumerate(hits.begin(), hits.end(), + + ccr::enumerate(hits.begin(), hits.end(), [&m, &rings, sd](HitResult &hit, size_t i) { - + // Point on the circle on the pin sphere Vec3d ps = rings.pinring(i); // This is the point on the circle on the back sphere Vec3d p = rings.backring(i); - + // Point ps is not on mesh but can be inside or // outside as well. This would cause many problems // with ray-casting. To detect the position we will // use the ray-casting result (which has an is_inside - // predicate). - + // predicate). + Vec3d n = (p - ps).normalized(); auto q = m.query_ray_hit(ps + sd * n, n); - + if (q.is_inside()) { // the hit is inside the model if (q.distance() > rings.rpin) { // If we are inside the model and the hit @@ -243,7 +250,7 @@ IndexedMesh::hit_result SupportTreeBuildsteps::pinhead_mesh_intersect( } else hit = q; }); - + return min_hit(hits); } @@ -252,20 +259,20 @@ IndexedMesh::hit_result SupportTreeBuildsteps::bridge_mesh_intersect( { static const size_t SAMPLES = 8; PointRing<SAMPLES> ring{dir}; - + using Hit = IndexedMesh::hit_result; - + // Hit results std::array<Hit, SAMPLES> hits; - - ccr::enumerate(hits.begin(), hits.end(), + + ccr::enumerate(hits.begin(), hits.end(), [this, r, src, /*ins_check,*/ &ring, dir, sd] (Hit &hit, size_t i) { // Point on the circle on the pin sphere Vec3d p = ring.get(i, src, r + sd); - + auto hr = m_mesh.query_ray_hit(p + r * dir, dir); - + if(/*ins_check && */hr.is_inside()) { if(hr.distance() > 2 * r + sd) hit = Hit(0.0); else { @@ -274,7 +281,7 @@ IndexedMesh::hit_result SupportTreeBuildsteps::bridge_mesh_intersect( } } else hit = hr; }); - + return min_hit(hits); } @@ -288,61 +295,61 @@ bool SupportTreeBuildsteps::interconnect(const Pillar &pillar, // shorter pillar is too short to start a new bridge but the taller // pillar could still be bridged with the shorter one. bool was_connected = false; - + Vec3d supper = pillar.startpoint(); Vec3d slower = nextpillar.startpoint(); Vec3d eupper = pillar.endpoint(); Vec3d elower = nextpillar.endpoint(); - + double zmin = m_builder.ground_level + m_cfg.base_height_mm; eupper(Z) = std::max(eupper(Z), zmin); elower(Z) = std::max(elower(Z), zmin); - + // The usable length of both pillars should be positive if(slower(Z) - elower(Z) < 0) return false; if(supper(Z) - eupper(Z) < 0) return false; - + double pillar_dist = distance(Vec2d{slower(X), slower(Y)}, Vec2d{supper(X), supper(Y)}); double bridge_distance = pillar_dist / std::cos(-m_cfg.bridge_slope); double zstep = pillar_dist * std::tan(-m_cfg.bridge_slope); - + if(pillar_dist < 2 * m_cfg.head_back_radius_mm || pillar_dist > m_cfg.max_pillar_link_distance_mm) return false; - + if(supper(Z) < slower(Z)) supper.swap(slower); if(eupper(Z) < elower(Z)) eupper.swap(elower); - + double startz = 0, endz = 0; - + startz = slower(Z) - zstep < supper(Z) ? slower(Z) - zstep : slower(Z); endz = eupper(Z) + zstep > elower(Z) ? eupper(Z) + zstep : eupper(Z); - + if(slower(Z) - eupper(Z) < std::abs(zstep)) { // no space for even one cross - + // Get max available space startz = std::min(supper(Z), slower(Z) - zstep); endz = std::max(eupper(Z) + zstep, elower(Z)); - + // Align to center double available_dist = (startz - endz); double rounds = std::floor(available_dist / std::abs(zstep)); startz -= 0.5 * (available_dist - rounds * std::abs(zstep)); } - + auto pcm = m_cfg.pillar_connection_mode; bool docrosses = pcm == PillarConnectionMode::cross || (pcm == PillarConnectionMode::dynamic && pillar_dist > 2*m_cfg.base_radius_mm); - + // 'sj' means starting junction, 'ej' is the end junction of a bridge. // They will be swapped in every iteration thus the zig-zag pattern. // According to a config parameter, a second bridge may be added which // results in a cross connection between the pillars. Vec3d sj = supper, ej = slower; sj(Z) = startz; ej(Z) = sj(Z) + zstep; - + // TODO: This is a workaround to not have a faulty last bridge while(ej(Z) >= eupper(Z) /*endz*/) { if(bridge_mesh_distance(sj, dirv(sj, ej), pillar.r) >= bridge_distance) @@ -350,7 +357,7 @@ bool SupportTreeBuildsteps::interconnect(const Pillar &pillar, m_builder.add_crossbridge(sj, ej, pillar.r); was_connected = true; } - + // double bridging: (crosses) if(docrosses) { Vec3d sjback(ej(X), ej(Y), sj(Z)); @@ -363,11 +370,11 @@ bool SupportTreeBuildsteps::interconnect(const Pillar &pillar, was_connected = true; } } - + sj.swap(ej); ej(Z) = sj(Z) + zstep; } - + return was_connected; } @@ -377,67 +384,67 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head, auto nearpillar = [this, nearpillar_id]() -> const Pillar& { return m_builder.pillar(nearpillar_id); }; - - if (m_builder.bridgecount(nearpillar()) > m_cfg.max_bridges_on_pillar) + + if (m_builder.bridgecount(nearpillar()) > m_cfg.max_bridges_on_pillar) return false; - + Vec3d headjp = head.junction_point(); Vec3d nearjp_u = nearpillar().startpoint(); Vec3d nearjp_l = nearpillar().endpoint(); - + double r = head.r_back_mm; double d2d = distance(to_2d(headjp), to_2d(nearjp_u)); double d3d = distance(headjp, nearjp_u); - + double hdiff = nearjp_u(Z) - headjp(Z); double slope = std::atan2(hdiff, d2d); - + Vec3d bridgestart = headjp; Vec3d bridgeend = nearjp_u; double max_len = r * m_cfg.max_bridge_length_mm / m_cfg.head_back_radius_mm; double max_slope = m_cfg.bridge_slope; double zdiff = 0.0; - + // check the default situation if feasible for a bridge if(d3d > max_len || slope > -max_slope) { // not feasible to connect the two head junctions. We have to search // for a suitable touch point. - + double Zdown = headjp(Z) + d2d * std::tan(-max_slope); Vec3d touchjp = bridgeend; touchjp(Z) = Zdown; double D = distance(headjp, touchjp); zdiff = Zdown - nearjp_u(Z); - + if(zdiff > 0) { Zdown -= zdiff; bridgestart(Z) -= zdiff; touchjp(Z) = Zdown; - + double t = bridge_mesh_distance(headjp, DOWN, r); - + // We can't insert a pillar under the source head to connect // with the nearby pillar's starting junction if(t < zdiff) return false; } - + if(Zdown <= nearjp_u(Z) && Zdown >= nearjp_l(Z) && D < max_len) bridgeend(Z) = Zdown; else return false; } - + // There will be a minimum distance from the ground where the // bridge is allowed to connect. This is an empiric value. double minz = m_builder.ground_level + 4 * head.r_back_mm; if(bridgeend(Z) < minz) return false; - + double t = bridge_mesh_distance(bridgestart, dirv(bridgestart, bridgeend), r); - + // Cannot insert the bridge. (further search might not worth the hassle) if(t < distance(bridgestart, bridgeend)) return false; - + std::lock_guard<ccr::BlockingMutex> lk(m_bridge_mutex); - + if (m_builder.bridgecount(nearpillar()) < m_cfg.max_bridges_on_pillar) { // A partial pillar is needed under the starting head. if(zdiff > 0) { @@ -447,31 +454,59 @@ bool SupportTreeBuildsteps::connect_to_nearpillar(const Head &head, } else { m_builder.add_bridge(head.id, bridgeend); } - + m_builder.increment_bridges(nearpillar()); } else return false; - + return true; } -bool SupportTreeBuildsteps::create_ground_pillar(const Vec3d &jp, +bool SupportTreeBuildsteps::create_ground_pillar(const Vec3d &hjp, const Vec3d &sourcedir, double radius, long head_id) { - double sd = m_cfg.pillar_base_safety_distance_mm; + Vec3d jp = hjp, endp = jp, dir = sourcedir; long pillar_id = SupportTreeNode::ID_UNSET; - bool can_add_base = radius >= m_cfg.head_back_radius_mm; - double base_r = can_add_base ? m_cfg.base_radius_mm : 0.; - double gndlvl = m_builder.ground_level; - if (!can_add_base) gndlvl -= m_mesh.ground_level_offset(); - Vec3d endp = {jp(X), jp(Y), gndlvl}; - double min_dist = sd + base_r + EPSILON; - bool normal_mode = true; - Vec3d dir = sourcedir; + bool can_add_base = false, non_head = false; + + double gndlvl = 0.; // The Z level where pedestals should be + double jp_gnd = 0.; // The lowest Z where a junction center can be + double gap_dist = 0.; // The gap distance between the model and the pad auto to_floor = [&gndlvl](const Vec3d &p) { return Vec3d{p.x(), p.y(), gndlvl}; }; + auto eval_limits = [this, &radius, &can_add_base, &gndlvl, &gap_dist, &jp_gnd] + (bool base_en = true) + { + can_add_base = base_en && radius >= m_cfg.head_back_radius_mm; + double base_r = can_add_base ? m_cfg.base_radius_mm : 0.; + gndlvl = m_builder.ground_level; + if (!can_add_base) gndlvl -= m_mesh.ground_level_offset(); + jp_gnd = gndlvl + (can_add_base ? 0. : m_cfg.head_back_radius_mm); + gap_dist = m_cfg.pillar_base_safety_distance_mm + base_r + EPSILON; + }; + + eval_limits(); + + // We are dealing with a mini pillar that's potentially too long + if (radius < m_cfg.head_back_radius_mm && jp.z() - gndlvl > 20 * radius) + { + std::optional<DiffBridge> diffbr = + search_widening_path(jp, dir, radius, m_cfg.head_back_radius_mm); + + if (diffbr && diffbr->endp.z() > jp_gnd) { + auto &br = m_builder.add_diffbridge(diffbr.value()); + if (head_id >= 0) m_builder.head(head_id).bridge_id = br.id; + endp = diffbr->endp; + radius = diffbr->end_r; + m_builder.add_junction(endp, radius); + non_head = true; + dir = diffbr->get_dir(); + eval_limits(); + } else return false; + } + if (m_cfg.object_elevation_mm < EPSILON) { // get a suitable direction for the corrector bridge. It is the @@ -479,101 +514,118 @@ bool SupportTreeBuildsteps::create_ground_pillar(const Vec3d &jp, // configured bridge slope. auto [polar, azimuth] = dir_to_spheric(dir); polar = PI - m_cfg.bridge_slope; - Vec3d dir = spheric_to_dir(polar, azimuth).normalized(); + Vec3d d = spheric_to_dir(polar, azimuth).normalized(); + double t = bridge_mesh_distance(endp, dir, radius); + double tmax = std::min(m_cfg.max_bridge_length_mm, t); + t = 0.; - // Check the distance of the endpoint and the closest point on model - // body. It should be greater than the min_dist which is - // the safety distance from the model. It includes the pad gap if in - // zero elevation mode. - // - // Try to move along the established bridge direction to dodge the - // forbidden region for the endpoint. - double t = -radius; - bool succ = true; - while (std::sqrt(m_mesh.squared_distance(to_floor(endp))) < min_dist || - !std::isinf(bridge_mesh_distance(endp, DOWN, radius))) { + double zd = endp.z() - jp_gnd; + double tmax2 = zd / std::sqrt(1 - m_cfg.bridge_slope * m_cfg.bridge_slope); + tmax = std::min(tmax, tmax2); + + Vec3d nexp = endp; + double dlast = 0.; + while (((dlast = std::sqrt(m_mesh.squared_distance(to_floor(nexp)))) < gap_dist || + !std::isinf(bridge_mesh_distance(nexp, DOWN, radius))) && t < tmax) { t += radius; - endp = jp + t * dir; - normal_mode = false; + nexp = endp + t * d; + } - if (t > m_cfg.max_bridge_length_mm || endp(Z) < gndlvl) { - if (head_id >= 0) m_builder.add_pillar(head_id, 0.); - succ = false; - break; + if (dlast < gap_dist && can_add_base) { + nexp = endp; + t = 0.; + can_add_base = false; + eval_limits(can_add_base); + + zd = endp.z() - jp_gnd; + tmax2 = zd / std::sqrt(1 - m_cfg.bridge_slope * m_cfg.bridge_slope); + tmax = std::min(tmax, tmax2); + + while (((dlast = std::sqrt(m_mesh.squared_distance(to_floor(nexp)))) < gap_dist || + !std::isinf(bridge_mesh_distance(nexp, DOWN, radius))) && t < tmax) { + t += radius; + nexp = endp + t * d; } } - if (!succ) { - if (can_add_base) { - can_add_base = false; - base_r = 0.; - gndlvl -= m_mesh.ground_level_offset(); - min_dist = sd + base_r + EPSILON; - endp = {jp(X), jp(Y), gndlvl + radius}; + // Could not find a path to avoid the pad gap + if (dlast < gap_dist) return false; - t = -radius; - while (std::sqrt(m_mesh.squared_distance(to_floor(endp))) < min_dist || - !std::isinf(bridge_mesh_distance(endp, DOWN, radius))) { - t += radius; - endp = jp + t * dir; - normal_mode = false; + if (t > 0.) { // Need to make additional bridge + const Bridge& br = m_builder.add_bridge(endp, nexp, radius); + if (head_id >= 0) m_builder.head(head_id).bridge_id = br.id; - if (t > m_cfg.max_bridge_length_mm || endp(Z) < (gndlvl + radius)) { - if (head_id >= 0) m_builder.add_pillar(head_id, 0.); - return false; - } - } - } else return false; + m_builder.add_junction(nexp, radius); + endp = nexp; + non_head = true; } } - double h = (jp - endp).norm(); + Vec3d gp = to_floor(endp); + double h = endp.z() - gp.z(); - // Check if the deduced route is sane and exit with error if not. - if (bridge_mesh_distance(jp, dir, radius) < h) { - if (head_id >= 0) m_builder.add_pillar(head_id, 0.); - return false; - } + pillar_id = head_id >= 0 && !non_head ? m_builder.add_pillar(head_id, h) : + m_builder.add_pillar(gp, h, radius); - // Straigh path down, no area to dodge - if (normal_mode) { - pillar_id = head_id >= 0 ? m_builder.add_pillar(head_id, h) : - m_builder.add_pillar(endp, h, radius); - - if (can_add_base) - add_pillar_base(pillar_id); - } else { - - // Insert the bridge to get around the forbidden area - Vec3d pgnd{endp.x(), endp.y(), gndlvl}; - pillar_id = m_builder.add_pillar(pgnd, endp.z() - gndlvl, radius); - - if (can_add_base) - add_pillar_base(pillar_id); - - m_builder.add_bridge(jp, endp, radius); - m_builder.add_junction(endp, radius); - - // Add a degenerated pillar and the bridge. - // The degenerate pillar will have zero length and it will - // prevent from queries of head_pillar() to have non-existing - // pillar when the head should have one. - if (head_id >= 0) - m_builder.add_pillar(head_id, 0.); - } + if (can_add_base) + add_pillar_base(pillar_id); if(pillar_id >= 0) // Save the pillar endpoint in the spatial index - m_pillar_index.guarded_insert(endp, unsigned(pillar_id)); + m_pillar_index.guarded_insert(m_builder.pillar(pillar_id).endpt, + unsigned(pillar_id)); return true; } +std::optional<DiffBridge> SupportTreeBuildsteps::search_widening_path( + const Vec3d &jp, const Vec3d &dir, double radius, double new_radius) +{ + double w = radius + 2 * m_cfg.head_back_radius_mm; + double stopval = w + jp.z() - m_builder.ground_level; + Optimizer<AlgNLoptSubplex> solver(get_criteria(m_cfg).stop_score(stopval)); + + auto [polar, azimuth] = dir_to_spheric(dir); + + double fallback_ratio = radius / m_cfg.head_back_radius_mm; + + auto oresult = solver.to_max().optimize( + [this, jp, radius, new_radius](double plr, double azm, double t) { + auto d = spheric_to_dir(plr, azm).normalized(); + double ret = pinhead_mesh_intersect(jp, d, radius, new_radius, t) + .distance(); + double down = bridge_mesh_distance(jp + t * d, d, new_radius); + + if (ret > t && std::isinf(down)) + ret += jp.z() - m_builder.ground_level; + + return ret; + }, + initvals({polar, azimuth, w}), // start with what we have + bounds({ + {PI - m_cfg.bridge_slope, PI}, // Must not exceed the slope limit + {-PI, PI}, // azimuth can be a full search + {radius + m_cfg.head_back_radius_mm, + fallback_ratio * m_cfg.max_bridge_length_mm} + })); + + if (oresult.score >= stopval) { + polar = std::get<0>(oresult.optimum); + azimuth = std::get<1>(oresult.optimum); + double t = std::get<2>(oresult.optimum); + Vec3d endp = jp + t * spheric_to_dir(polar, azimuth); + + return DiffBridge(jp, endp, radius, m_cfg.head_back_radius_mm); + } + + return {}; +} + void SupportTreeBuildsteps::filter() { // Get the points that are too close to each other and keep only the // first one auto aliases = cluster(m_points, D_SP, 2); - + PtIndices filtered_indices; filtered_indices.reserve(aliases.size()); m_iheads.reserve(aliases.size()); @@ -582,49 +634,62 @@ void SupportTreeBuildsteps::filter() // Here we keep only the front point of the cluster. filtered_indices.emplace_back(a.front()); } - + // calculate the normals to the triangles for filtered points auto nmls = sla::normals(m_points, m_mesh, m_cfg.head_front_radius_mm, m_thr, filtered_indices); - + // Not all of the support points have to be a valid position for // support creation. The angle may be inappropriate or there may // not be enough space for the pinhead. Filtering is applied for // these reasons. - - ccr::SpinningMutex mutex; - auto addfn = [&mutex](PtIndices &container, unsigned val) { - std::lock_guard<ccr::SpinningMutex> lk(mutex); - container.emplace_back(val); - }; - - auto filterfn = [this, &nmls, addfn](unsigned fidx, size_t i) { + + std::vector<Head> heads; heads.reserve(m_support_pts.size()); + for (const SupportPoint &sp : m_support_pts) { m_thr(); - + heads.emplace_back( + std::nan(""), + sp.head_front_radius, + 0., + m_cfg.head_penetration_mm, + Vec3d::Zero(), // dir + sp.pos.cast<double>() // displacement + ); + } + + std::function<void(unsigned, size_t, double)> filterfn; + filterfn = [this, &nmls, &heads, &filterfn](unsigned fidx, size_t i, double back_r) { + m_thr(); + auto n = nmls.row(Eigen::Index(i)); - + // for all normals we generate the spherical coordinates and // saturate the polar angle to 45 degrees from the bottom then // convert back to standard coordinates to get the new normal. // Then we just create a quaternion from the two normals // (Quaternion::FromTwoVectors) and apply the rotation to the // arrow head. - + auto [polar, azimuth] = dir_to_spheric(n); - + // skip if the tilt is not sane - if(polar < PI - m_cfg.normal_cutoff_angle) return; - + if (polar < PI - m_cfg.normal_cutoff_angle) return; + // We saturate the polar angle to 3pi/4 - polar = std::max(polar, 3*PI / 4); + polar = std::max(polar, PI - m_cfg.bridge_slope); // save the head (pinpoint) position Vec3d hp = m_points.row(fidx); + double lmin = m_cfg.head_width_mm, lmax = lmin; + + if (back_r < m_cfg.head_back_radius_mm) { + lmin = 0., lmax = m_cfg.head_penetration_mm; + } + // The distance needed for a pinhead to not collide with model. - double w = m_cfg.head_width_mm + - m_cfg.head_back_radius_mm + - 2*m_cfg.head_front_radius_mm; + double w = lmin + 2 * back_r + 2 * m_cfg.head_front_radius_mm - + m_cfg.head_penetration_mm; double pin_r = double(m_support_pts[fidx].head_front_radius); @@ -632,113 +697,69 @@ void SupportTreeBuildsteps::filter() auto nn = spheric_to_dir(polar, azimuth).normalized(); // check available distance - IndexedMesh::hit_result t - = pinhead_mesh_intersect(hp, // touching point - nn, // normal - pin_r, - m_cfg.head_back_radius_mm, - w); - - if(t.distance() <= w) { + IndexedMesh::hit_result t = pinhead_mesh_intersect(hp, nn, pin_r, + back_r, w); + if (t.distance() < w) { // Let's try to optimize this angle, there might be a // viable normal that doesn't collide with the model // geometry and its very close to the default. - StopCriteria stc; - stc.max_iterations = m_cfg.optimizer_max_iterations; - stc.relative_score_difference = m_cfg.optimizer_rel_score_diff; - stc.stop_score = w; // space greater than w is enough - GeneticOptimizer solver(stc); - solver.seed(0); // we want deterministic behavior + // stc.stop_score = w; // space greater than w is enough + Optimizer<AlgNLoptGenetic> solver(get_criteria(m_cfg)); + solver.seed(0); + //solver.seed(0); // we want deterministic behavior - auto oresult = solver.optimize_max( - [this, pin_r, w, hp](double plr, double azm) + auto oresult = solver.to_max().optimize( + [this, pin_r, back_r, hp](double plr, double azm, double l) { auto dir = spheric_to_dir(plr, azm).normalized(); double score = pinhead_mesh_intersect( - hp, dir, pin_r, m_cfg.head_back_radius_mm, w).distance(); + hp, dir, pin_r, back_r, l).distance(); return score; }, - initvals(polar, azimuth), // start with what we have - bound(3 * PI / 4, PI), // Must not exceed the tilt limit - bound(-PI, PI) // azimuth can be a full search - ); + initvals({polar, azimuth, (lmin + lmax) / 2.}), // start with what we have + bounds({ + {PI - m_cfg.bridge_slope, PI}, // Must not exceed the tilt limit + {-PI, PI}, // azimuth can be a full search + {lmin, lmax} + })); if(oresult.score > w) { polar = std::get<0>(oresult.optimum); azimuth = std::get<1>(oresult.optimum); nn = spheric_to_dir(polar, azimuth).normalized(); + lmin = std::get<2>(oresult.optimum); t = IndexedMesh::hit_result(oresult.score); } } - // save the verified and corrected normal - m_support_nmls.row(fidx) = nn; + if (t.distance() > w && hp(Z) + w * nn(Z) >= m_builder.ground_level) { + Head &h = heads[fidx]; + h.id = fidx; h.dir = nn; h.width_mm = lmin; h.r_back_mm = back_r; + } else if (back_r > m_cfg.head_fallback_radius_mm) { + filterfn(fidx, i, m_cfg.head_fallback_radius_mm); + } + }; - if (t.distance() > w) { - // Check distance from ground, we might have zero elevation. - if (hp(Z) + w * nn(Z) < m_builder.ground_level) { - addfn(m_iheadless, fidx); - } else { - // mark the point for needing a head. - addfn(m_iheads, fidx); - } - } else if (polar >= 3 * PI / 4) { - // Headless supports do not tilt like the headed ones - // so the normal should point almost to the ground. - addfn(m_iheadless, fidx); + ccr::enumerate(filtered_indices.begin(), filtered_indices.end(), + [this, &filterfn](unsigned fidx, size_t i) { + filterfn(fidx, i, m_cfg.head_back_radius_mm); + }); + + for (size_t i = 0; i < heads.size(); ++i) + if (heads[i].is_valid()) { + m_builder.add_head(i, heads[i]); + m_iheads.emplace_back(i); } - }; - - ccr::enumerate(filtered_indices.begin(), filtered_indices.end(), filterfn); - m_thr(); } void SupportTreeBuildsteps::add_pinheads() { - for (unsigned i : m_iheads) { - m_thr(); - m_builder.add_head( - i, - m_cfg.head_back_radius_mm, - m_support_pts[i].head_front_radius, - m_cfg.head_width_mm, - m_cfg.head_penetration_mm, - m_support_nmls.row(i), // dir - m_support_pts[i].pos.cast<double>() // displacement - ); - } - - for (unsigned i : m_iheadless) { - const auto R = double(m_support_pts[i].head_front_radius); - - // The support point position on the mesh - Vec3d sph = m_support_pts[i].pos.cast<double>(); - - // Get an initial normal from the filtering step - Vec3d n = m_support_nmls.row(i); - - // First we need to determine the available space for a mini pinhead. - // The goal is the move away from the model a little bit to make the - // contact point small as possible and avoid pearcing the model body. - double back_r = m_cfg.head_fallback_radius_mm; - double max_w = 2 * R; - double pin_space = std::min(max_w, - pinhead_mesh_intersect(sph, n, R, back_r, - max_w, 0.) - .distance()); - - if (pin_space <= 0) continue; - - m_iheads.emplace_back(i); - m_builder.add_head(i, back_r, R, pin_space, - m_cfg.head_penetration_mm, n, sph); - } } void SupportTreeBuildsteps::classify() @@ -747,37 +768,37 @@ void SupportTreeBuildsteps::classify() PtIndices ground_head_indices; ground_head_indices.reserve(m_iheads.size()); m_iheads_onmodel.reserve(m_iheads.size()); - + // First we decide which heads reach the ground and can be full // pillars and which shall be connected to the model surface (or // search a suitable path around the surface that leads to the // ground -- TODO) for(unsigned i : m_iheads) { m_thr(); - - auto& head = m_builder.head(i); + + Head &head = m_builder.head(i); double r = head.r_back_mm; Vec3d headjp = head.junction_point(); - + // collision check auto hit = bridge_mesh_intersect(headjp, DOWN, r); - + if(std::isinf(hit.distance())) ground_head_indices.emplace_back(i); else if(m_cfg.ground_facing_only) head.invalidate(); else m_iheads_onmodel.emplace_back(i); - + m_head_to_ground_scans[i] = hit; } - + // We want to search for clusters of points that are far enough // from each other in the XY plane to not cross their pillar bases // These clusters of support points will join in one pillar, // possibly in their centroid support point. - + auto pointfn = [this](unsigned i) { return m_builder.head(i).junction_point(); }; - + auto predicate = [this](const PointIndexEl &e1, const PointIndexEl &e2) { double d2d = distance(to_2d(e1.first), to_2d(e2.first)); @@ -794,10 +815,10 @@ void SupportTreeBuildsteps::routing_to_ground() { ClusterEl cl_centroids; cl_centroids.reserve(m_pillar_clusters.size()); - + for (auto &cl : m_pillar_clusters) { m_thr(); - + // place all the centroid head positions into the index. We // will query for alternative pillar positions. If a sidehead // cannot connect to the cluster centroid, we have to search @@ -805,9 +826,9 @@ void SupportTreeBuildsteps::routing_to_ground() // elements in the cluster, the centroid is arbitrary and the // sidehead is allowed to connect to a nearby pillar to // increase structural stability. - + if (cl.empty()) continue; - + // get the current cluster centroid auto & thr = m_thr; const auto &points = m_points; @@ -821,11 +842,11 @@ void SupportTreeBuildsteps::routing_to_ground() assert(lcid >= 0); unsigned hid = cl[size_t(lcid)]; // Head ID - + cl_centroids.emplace_back(hid); - + Head &h = m_builder.head(hid); - + if (!create_ground_pillar(h.junction_point(), h.dir, h.r_back_mm, h.id)) { BOOST_LOG_TRIVIAL(warning) << "Pillar cannot be created for support point id: " << hid; @@ -833,34 +854,32 @@ void SupportTreeBuildsteps::routing_to_ground() continue; } } - + // now we will go through the clusters ones again and connect the // sidepoints with the cluster centroid (which is a ground pillar) // or a nearby pillar if the centroid is unreachable. size_t ci = 0; for (auto cl : m_pillar_clusters) { m_thr(); - + auto cidx = cl_centroids[ci++]; - - // TODO: don't consider the cluster centroid but calculate a - // central position where the pillar can be placed. this way - // the weight is distributed more effectively on the pillar. - - auto centerpillarID = m_builder.head_pillar(cidx).id; - - for (auto c : cl) { - m_thr(); - if (c == cidx) continue; - - auto &sidehead = m_builder.head(c); - - if (!connect_to_nearpillar(sidehead, centerpillarID) && - !search_pillar_and_connect(sidehead)) { - Vec3d pstart = sidehead.junction_point(); - // Vec3d pend = Vec3d{pstart(X), pstart(Y), gndlvl}; - // Could not find a pillar, create one - create_ground_pillar(pstart, sidehead.dir, sidehead.r_back_mm, sidehead.id); + + auto q = m_pillar_index.query(m_builder.head(cidx).junction_point(), 1); + if (!q.empty()) { + long centerpillarID = q.front().second; + for (auto c : cl) { + m_thr(); + if (c == cidx) continue; + + auto &sidehead = m_builder.head(c); + + if (!connect_to_nearpillar(sidehead, centerpillarID) && + !search_pillar_and_connect(sidehead)) { + Vec3d pstart = sidehead.junction_point(); + // Vec3d pend = Vec3d{pstart(X), pstart(Y), gndlvl}; + // Could not find a pillar, create one + create_ground_pillar(pstart, sidehead.dir, sidehead.r_back_mm, sidehead.id); + } } } } @@ -876,9 +895,9 @@ bool SupportTreeBuildsteps::connect_to_ground(Head &head, const Vec3d &dir) while (d < t && !std::isinf(tdown = bridge_mesh_distance(hjp + d * dir, DOWN, r))) d += r; - + if(!std::isinf(tdown)) return false; - + Vec3d endp = hjp + d * dir; bool ret = false; @@ -886,38 +905,33 @@ bool SupportTreeBuildsteps::connect_to_ground(Head &head, const Vec3d &dir) m_builder.add_bridge(head.id, endp); m_builder.add_junction(endp, head.r_back_mm); } - + return ret; } bool SupportTreeBuildsteps::connect_to_ground(Head &head) { if (connect_to_ground(head, head.dir)) return true; - + // Optimize bridge direction: // Straight path failed so we will try to search for a suitable // direction out of the cavity. auto [polar, azimuth] = dir_to_spheric(head.dir); - - StopCriteria stc; - stc.max_iterations = m_cfg.optimizer_max_iterations; - stc.relative_score_difference = m_cfg.optimizer_rel_score_diff; - stc.stop_score = 1e6; - GeneticOptimizer solver(stc); + + Optimizer<AlgNLoptGenetic> solver(get_criteria(m_cfg).stop_score(1e6)); solver.seed(0); // we want deterministic behavior - + double r_back = head.r_back_mm; - Vec3d hjp = head.junction_point(); - auto oresult = solver.optimize_max( + Vec3d hjp = head.junction_point(); + auto oresult = solver.to_max().optimize( [this, hjp, r_back](double plr, double azm) { Vec3d n = spheric_to_dir(plr, azm).normalized(); return bridge_mesh_distance(hjp, n, r_back); }, - initvals(polar, azimuth), // let's start with what we have - bound(3*PI/4, PI), // Must not exceed the slope limit - bound(-PI, PI) // azimuth can be a full range search - ); - + initvals({polar, azimuth}), // let's start with what we have + bounds({ {PI - m_cfg.bridge_slope, PI}, {-PI, PI} }) + ); + Vec3d bridgedir = spheric_to_dir(oresult.optimum).normalized(); return connect_to_ground(head, bridgedir); } @@ -925,10 +939,10 @@ bool SupportTreeBuildsteps::connect_to_ground(Head &head) bool SupportTreeBuildsteps::connect_to_model_body(Head &head) { if (head.id <= SupportTreeNode::ID_UNSET) return false; - + auto it = m_head_to_ground_scans.find(unsigned(head.id)); if (it == m_head_to_ground_scans.end()) return false; - + auto &hit = it->second; if (!hit.is_hit()) { @@ -943,9 +957,11 @@ bool SupportTreeBuildsteps::connect_to_model_body(Head &head) // The width of the tail head that we would like to have... h = std::min(hit.distance() - head.r_back_mm, h); - - if(h <= 0.) return false; - + + // If this is a mini pillar dont bother with the tail width, can be 0. + if (head.r_back_mm < m_cfg.head_back_radius_mm) h = std::max(h, 0.); + else if (h <= 0.) return false; + Vec3d endp{hjp(X), hjp(Y), hjp(Z) - hit.distance() + h}; auto center_hit = m_mesh.query_ray_hit(hjp, DOWN); @@ -969,7 +985,7 @@ bool SupportTreeBuildsteps::connect_to_model_body(Head &head) m_cfg.head_penetration_mm, taildir, hitp); m_pillar_index.guarded_insert(pill.endpoint(), pill.id); - + return true; } @@ -1011,7 +1027,7 @@ bool SupportTreeBuildsteps::search_pillar_and_connect(const Head &source) } void SupportTreeBuildsteps::routing_to_model() -{ +{ // We need to check if there is an easy way out to the bed surface. // If it can be routed there with a bridge shorter than // min_bridge_distance. @@ -1019,23 +1035,23 @@ void SupportTreeBuildsteps::routing_to_model() ccr::enumerate(m_iheads_onmodel.begin(), m_iheads_onmodel.end(), [this] (const unsigned idx, size_t) { m_thr(); - + auto& head = m_builder.head(idx); - + // Search nearby pillar if (search_pillar_and_connect(head)) { return; } - + // Cannot connect to nearby pillar. We will try to search for // a route to the ground. if (connect_to_ground(head)) { return; } - + // No route to the ground, so connect to the model body as a last resort if (connect_to_model_body(head)) { return; } - + // We have failed to route this head. BOOST_LOG_TRIVIAL(warning) << "Failed to route model facing support point. ID: " << idx; - + head.invalidate(); }); } @@ -1045,19 +1061,19 @@ void SupportTreeBuildsteps::interconnect_pillars() // Now comes the algorithm that connects pillars with each other. // Ideally every pillar should be connected with at least one of its // neighbors if that neighbor is within max_pillar_link_distance - + // Pillars with height exceeding H1 will require at least one neighbor // to connect with. Height exceeding H2 require two neighbors. double H1 = m_cfg.max_solo_pillar_height_mm; double H2 = m_cfg.max_dual_pillar_height_mm; double d = m_cfg.max_pillar_link_distance_mm; - + //A connection between two pillars only counts if the height ratio is // bigger than 50% double min_height_ratio = 0.5; - + std::set<unsigned long> pairs; - + // A function to connect one pillar with its neighbors. THe number of // neighbors is given in the configuration. This function if called // for every pillar in the pillar index. A pair of pillar will not @@ -1067,68 +1083,68 @@ void SupportTreeBuildsteps::interconnect_pillars() [this, d, &pairs, min_height_ratio, H1] (const PointIndexEl& el) { Vec3d qp = el.first; // endpoint of the pillar - + const Pillar& pillar = m_builder.pillar(el.second); // actual pillar - + // Get the max number of neighbors a pillar should connect to unsigned neighbors = m_cfg.pillar_cascade_neighbors; - + // connections are already enough for the pillar if(pillar.links >= neighbors) return; - + double max_d = d * pillar.r / m_cfg.head_back_radius_mm; // Query all remaining points within reach auto qres = m_pillar_index.query([qp, max_d](const PointIndexEl& e){ return distance(e.first, qp) < max_d; }); - + // sort the result by distance (have to check if this is needed) std::sort(qres.begin(), qres.end(), [qp](const PointIndexEl& e1, const PointIndexEl& e2){ return distance(e1.first, qp) < distance(e2.first, qp); }); - + for(auto& re : qres) { // process the queried neighbors - + if(re.second == el.second) continue; // Skip self - + auto a = el.second, b = re.second; - + // Get unique hash for the given pair (order doesn't matter) auto hashval = pairhash(a, b); - + // Search for the pair amongst the remembered pairs if(pairs.find(hashval) != pairs.end()) continue; - + const Pillar& neighborpillar = m_builder.pillar(re.second); - + // this neighbor is occupied, skip if (neighborpillar.links >= neighbors) continue; if (neighborpillar.r < pillar.r) continue; - + if(interconnect(pillar, neighborpillar)) { pairs.insert(hashval); - + // If the interconnection length between the two pillars is // less than 50% of the longer pillar's height, don't count if(pillar.height < H1 || neighborpillar.height / pillar.height > min_height_ratio) m_builder.increment_links(pillar); - + if(neighborpillar.height < H1 || pillar.height / neighborpillar.height > min_height_ratio) m_builder.increment_links(neighborpillar); - + } - + // connections are enough for one pillar if(pillar.links >= neighbors) break; } }; - + // Run the cascade for the pillars in the index m_pillar_index.foreach(cascadefn); - + // We would be done here if we could allow some pillars to not be // connected with any neighbors. But this might leave the support tree // unprintable. @@ -1136,16 +1152,16 @@ void SupportTreeBuildsteps::interconnect_pillars() // The current solution is to insert additional pillars next to these // lonely pillars. One or even two additional pillar might get inserted // depending on the length of the lonely pillar. - + size_t pillarcount = m_builder.pillarcount(); - + // Again, go through all pillars, this time in the whole support tree // not just the index. for(size_t pid = 0; pid < pillarcount; pid++) { auto pillar = [this, pid]() { return m_builder.pillar(pid); }; - + // Decide how many additional pillars will be needed: - + unsigned needpillars = 0; if (pillar().bridges > m_cfg.max_bridges_on_pillar) needpillars = 3; @@ -1156,28 +1172,28 @@ void SupportTreeBuildsteps::interconnect_pillars() // No neighbors could be found and the pillar is too long. needpillars = 1; } - + needpillars = std::max(pillar().links, needpillars) - pillar().links; if (needpillars == 0) continue; - + // Search for new pillar locations: - + bool found = false; double alpha = 0; // goes to 2Pi double r = 2 * m_cfg.base_radius_mm; Vec3d pillarsp = pillar().startpoint(); - + // temp value for starting point detection Vec3d sp(pillarsp(X), pillarsp(Y), pillarsp(Z) - r); - + // A vector of bool for placement feasbility std::vector<bool> canplace(needpillars, false); std::vector<Vec3d> spts(needpillars); // vector of starting points - + double gnd = m_builder.ground_level; double min_dist = m_cfg.pillar_base_safety_distance_mm + m_cfg.base_radius_mm + EPSILON; - + while(!found && alpha < 2*PI) { for (unsigned n = 0; n < needpillars && (!n || canplace[n - 1]); @@ -1188,25 +1204,25 @@ void SupportTreeBuildsteps::interconnect_pillars() s(X) += std::cos(a) * r; s(Y) += std::sin(a) * r; spts[n] = s; - + // Check the path vertically down Vec3d check_from = s + Vec3d{0., 0., pillar().r}; auto hr = bridge_mesh_intersect(check_from, DOWN, pillar().r); Vec3d gndsp{s(X), s(Y), gnd}; - + // If the path is clear, check for pillar base collisions canplace[n] = std::isinf(hr.distance()) && std::sqrt(m_mesh.squared_distance(gndsp)) > min_dist; } - + found = std::all_of(canplace.begin(), canplace.end(), [](bool v) { return v; }); - + // 20 angles will be tried... alpha += 0.1 * PI; } - + std::vector<long> newpills; newpills.reserve(needpillars); @@ -1247,7 +1263,7 @@ void SupportTreeBuildsteps::interconnect_pillars() m_builder.increment_links(nxpll); } } - + m_pillar_index.foreach(cascadefn); } } diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.hpp b/src/libslic3r/SLA/SupportTreeBuildsteps.hpp index d19194a87..013666f07 100644 --- a/src/libslic3r/SLA/SupportTreeBuildsteps.hpp +++ b/src/libslic3r/SLA/SupportTreeBuildsteps.hpp @@ -45,6 +45,11 @@ inline Vec3d spheric_to_dir(const std::pair<double, double> &v) return spheric_to_dir(v.first, v.second); } +inline Vec3d spheric_to_dir(const std::array<double, 2> &v) +{ + return spheric_to_dir(v[0], v[1]); +} + // Give points on a 3D ring with given center, radius and orientation // method based on: // https://math.stackexchange.com/questions/73237/parametric-equation-of-a-circle-in-3d-space @@ -249,7 +254,8 @@ class SupportTreeBuildsteps { double width) { return pinhead_mesh_intersect(s, dir, r_pin, r_back, width, - m_cfg.safety_distance_mm); + r_back * m_cfg.safety_distance_mm / + m_cfg.head_back_radius_mm); } // Checking bridge (pillar and stick as well) intersection with the model. @@ -271,7 +277,9 @@ class SupportTreeBuildsteps { const Vec3d& dir, double r) { - return bridge_mesh_intersect(s, dir, r, m_cfg.safety_distance_mm); + return bridge_mesh_intersect(s, dir, r, + r * m_cfg.safety_distance_mm / + m_cfg.head_back_radius_mm); } template<class...Args> @@ -311,6 +319,11 @@ class SupportTreeBuildsteps { m_builder.add_pillar_base(pid, m_cfg.base_height_mm, m_cfg.base_radius_mm); } + std::optional<DiffBridge> search_widening_path(const Vec3d &jp, + const Vec3d &dir, + double radius, + double new_radius); + public: SupportTreeBuildsteps(SupportTreeBuilder & builder, const SupportableMesh &sm); diff --git a/src/libslic3r/SLA/SupportTreeMesher.hpp b/src/libslic3r/SLA/SupportTreeMesher.hpp index a086680c3..63182745d 100644 --- a/src/libslic3r/SLA/SupportTreeMesher.hpp +++ b/src/libslic3r/SLA/SupportTreeMesher.hpp @@ -94,6 +94,24 @@ inline Contour3D get_mesh(const Bridge &br, size_t steps) return mesh; } +inline Contour3D get_mesh(const DiffBridge &br, size_t steps) +{ + double h = br.get_length(); + Contour3D mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps); + + using Quaternion = Eigen::Quaternion<double>; + + // We rotate the head to the specified direction. The head's pointing + // side is facing upwards so this means that it would hold a support + // point with a normal pointing straight down. This is the reason of + // the -1 z coordinate + auto quatern = Quaternion::FromTwoVectors(Vec3d{0, 0, 1}, br.get_dir()); + + for(auto& p : mesh.points) p = quatern * p + br.startp; + + return mesh; +} + }} // namespace Slic3r::sla #endif // SUPPORTTREEMESHER_HPP diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp index 9a9c762e3..dad2b9097 100644 --- a/tests/sla_print/sla_print_tests.cpp +++ b/tests/sla_print/sla_print_tests.cpp @@ -4,6 +4,8 @@ #include "sla_test_utils.hpp" +#include <libslic3r/SLA/SupportTreeMesher.hpp> + namespace { const char *const BELOW_PAD_TEST_OBJECTS[] = { @@ -228,3 +230,12 @@ TEST_CASE("Triangle mesh conversions should be correct", "[SLAConversions]") cntr.from_obj(infile); } } + +TEST_CASE("halfcone test", "[halfcone]") { + sla::DiffBridge br{Vec3d{1., 1., 1.}, Vec3d{10., 10., 10.}, 0.25, 0.5}; + + TriangleMesh m = sla::to_triangle_mesh(sla::get_mesh(br, 45)); + + m.require_shared_vertices(); + m.WriteOBJFile("Halfcone.obj"); +}