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

@@ -203,6 +203,7 @@ add_library(libslic3r STATIC
     SimplifyMeshImpl.hpp
     SimplifyMesh.cpp
     MarchingSquares.hpp
+    Optimizer.hpp
     ${OpenVDBUtils_SOURCES}
     SLA/Pad.hpp
     SLA/Pad.cpp

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

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());

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));

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)
     {

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);
 

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

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");
+}