diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 09f75c747..e30811133 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -215,7 +215,9 @@ add_library(libslic3r STATIC
     SimplifyMeshImpl.hpp
     SimplifyMesh.cpp
     MarchingSquares.hpp
-    Optimizer.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
new file mode 100644
index 000000000..2daef538e
--- /dev/null
+++ b/src/libslic3r/Optimize/BruteforceOptimizer.hpp
@@ -0,0 +1,140 @@
+#ifndef BRUTEFORCEOPTIMIZER_HPP
+#define BRUTEFORCEOPTIMIZER_HPP
+
+#include <libslic3r/Optimize/Optimizer.hpp>
+
+namespace Slic3r { namespace opt {
+
+namespace detail {
+// Implementing a bruteforce optimizer
+
+// Return the number of iterations needed to reach a specific grid position (idx)
+template<size_t N>
+long num_iter(const std::array<size_t, N> &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;
+    size_t gridsz;
+
+    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<int D, size_t N, class Fn, class Cmp>
+    bool run(std::array<size_t, N> &idx,
+             Result<N> &result,
+             const Bounds<N> &bounds,
+             Fn &&fn,
+             Cmp &&cmp)
+    {
+        if (stc.stop_condition()) return false;
+
+        if constexpr (D < 0) { // Let's evaluate fn
+            Input<N> inp;
+
+            auto max_iter = stc.max_iterations();
+            if (max_iter && num_iter(idx, gridsz) >= max_iter)
+                return false;
+
+            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)) { // 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;
+            }
+
+        } else {
+            for (size_t i = 0; i < gridsz; ++i) {
+                idx[D] = i; // Mark the current grid position and dig down
+                if (!run<D - 1>(idx, result, bounds, std::forward<Fn>(fn),
+                                std::forward<Cmp>(cmp)))
+                    return false;
+            }
+        }
+
+        return true;
+    }
+
+    template<class Fn, size_t N>
+    Result<N> optimize(Fn&& fn,
+                       const Input<N> &/*initvals*/,
+                       const Bounds<N>& bounds)
+    {
+        std::array<size_t, N> idx = {};
+        Result<N> result;
+
+        if (to_min) {
+            result.score = std::numeric_limits<double>::max();
+            run<int(N) - 1>(idx, result, bounds, std::forward<Fn>(fn),
+                            std::less<double>{});
+        }
+        else {
+            result.score = std::numeric_limits<double>::lowest();
+            run<int(N) - 1>(idx, result, bounds, std::forward<Fn>(fn),
+                            std::greater<double>{});
+        }
+
+        return result;
+    }
+};
+
+} // namespace detail
+
+using AlgBruteForce = detail::AlgBurteForce;
+
+template<>
+class Optimizer<AlgBruteForce> {
+    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<class Func, size_t N>
+    Result<N> optimize(Func&& func,
+                       const Input<N> &initvals,
+                       const Bounds<N>& bounds)
+    {
+        return m_alg.optimize(std::forward<Func>(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 6495ae7ff..826b1632a 100644
--- a/src/libslic3r/Optimizer.hpp
+++ b/src/libslic3r/Optimize/NLoptOptimizer.hpp
@@ -12,134 +12,11 @@
 #endif
 
 #include <utility>
-#include <tuple>
-#include <array>
-#include <cmath>
-#include <functional>
-#include <limits>
-#include <cassert>
+
+#include <libslic3r/Optimize/Optimizer.hpp>
 
 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 class to use optimization methods involving gradient.
-template<size_t N> struct ScoreGradient {
-    double score;
-    std::optional<std::array<double, N>> gradient;
-
-    ScoreGradient(double s, const std::array<double, N> &grad)
-        : score{s}, gradient{grad}
-    {}
-};
-
-// 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> &to_min() { return *this; }
-    Optimizer<Method> &to_max() { return *this; }
-    Optimizer<Method> &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
@@ -166,19 +43,6 @@ struct IsNLoptAlg<NLoptAlgComb<a1, a2>> {
 template<class M, class T = void>
 using NLoptOnly = std::enable_if_t<IsNLoptAlg<M>::value, T>;
 
-// Helper to convert C style array to std::array. The copy should be optimized
-// away with modern compilers.
-template<size_t N, class T> auto to_arr(const T *a)
-{
-    std::array<T, N> r;
-    std::copy(a, a + N, std::begin(r));
-    return r;
-}
-
-template<size_t N, class T> auto to_arr(const T (&a) [N])
-{
-    return to_arr<N>(static_cast<const T *>(a));
-}
 
 enum class OptDir { MIN, MAX }; // Where to optimize
 
@@ -357,23 +221,12 @@ public:
     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); }
-template<size_t N> auto score_gradient(double s, const double (&grad)[N])
-{
-    return ScoreGradient<N>(s, detail::to_arr(grad));
-}
-
-// Predefinded NLopt algorithms that are used in the codebase
+// Predefinded NLopt algorithms
 using AlgNLoptGenetic = detail::NLoptAlgComb<NLOPT_GN_ESCH>;
 using AlgNLoptSubplex = detail::NLoptAlg<NLOPT_LN_SBPLX>;
 using AlgNLoptSimplex = detail::NLoptAlg<NLOPT_LN_NELDERMEAD>;
-
-// TODO: define others if needed...
-
-// Helper defs for pre-crafted global and local optimizers that work well.
-using DefaultGlobalOptimizer = Optimizer<AlgNLoptGenetic>;
-using DefaultLocalOptimizer  = Optimizer<AlgNLoptSubplex>;
+using AlgNLoptDIRECT  = detail::NLoptAlg<NLOPT_GN_DIRECT>;
+using AlgNLoptMLSL    = detail::NLoptAlg<NLOPT_GN_MLSL>;
 
 }} // 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 <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;     // Method dependent
+    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 class to use optimization methods involving gradient.
+template<size_t N> struct ScoreGradient {
+    double score;
+    std::optional<std::array<double, N>> gradient;
+
+    ScoreGradient(double s, const std::array<double, N> &grad)
+        : score{s}, gradient{grad}
+    {}
+};
+
+// 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!");
+    }
+
+    // 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<N> &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<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 to convert C style array to std::array. The copy should be optimized
+// away with modern compilers.
+template<size_t N, class T> auto to_arr(const T *a)
+{
+    std::array<T, N> r;
+    std::copy(a, a + N, std::begin(r));
+    return r;
+}
+
+template<size_t N, class T> auto to_arr(const T (&a) [N])
+{
+    return to_arr<N>(static_cast<const T *>(a));
+}
+
+} // namespace detail
+
+// Helper functions to create bounds, initial value
+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); }
+template<size_t N> auto score_gradient(double s, const double (&grad)[N])
+{
+    return ScoreGradient<N>(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..300024c76 100644
--- a/src/libslic3r/SLA/Concurrency.hpp
+++ b/src/libslic3r/SLA/Concurrency.hpp
@@ -4,7 +4,11 @@
 #include <tbb/spin_mutex.h>
 #include <tbb/mutex.h>
 #include <tbb/parallel_for.h>
+#include <tbb/parallel_reduce.h>
+
 #include <algorithm>
+#include <numeric>
+
 #include <libslic3r/libslic3r.h>
 
 namespace Slic3r {
@@ -21,28 +25,56 @@ template<> struct _ccr<true>
     using SpinningMutex = tbb::spin_mutex;
     using BlockingMutex = tbb::mutex;
 
+    template<class Fn, class It>
+    static IteratorOnly<It, void> loop_(const tbb::blocked_range<It> &range, Fn &&fn)
+    {
+        for (auto &el : range) fn(el);
+    }
+
+    template<class Fn, class I>
+    static IntegerOnly<I, void> loop_(const tbb::blocked_range<I> &range, Fn &&fn)
+    {
+        for (I i = range.begin(); i < range.end(); ++i) fn(i);
+    }
+
     template<class It, class Fn>
-    static IteratorOnly<It, void> 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>(fn));
         });
     }
 
-    template<class I, class Fn>
-    static IntegerOnly<I, void> for_each(I      from,
-                                         I      to,
-                                         Fn &&  fn,
-                                         size_t granularity = 1)
+    template<class I, class MergeFn, class T, class AccessFn>
+    static T reduce(I          from,
+                    I          to,
+                    const T   &init,
+                    MergeFn  &&mergefn,
+                    AccessFn &&access,
+                    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, access(i)); });
+                return acc;
+            },
+            std::forward<MergeFn>(mergefn));
+    }
+
+    template<class I, class MergeFn, class T>
+    static IteratorOnly<I, T> reduce(I         from,
+                                     I         to,
+                                     const T & init,
+                                     MergeFn &&mergefn,
+                                     size_t    granularity = 1)
+    {
+        return reduce(
+            from, to, init, std::forward<MergeFn>(mergefn),
+            [](typename I::value_type &i) { return i; }, granularity);
     }
 };
 
@@ -55,23 +87,52 @@ public:
     using SpinningMutex = _Mtx;
     using BlockingMutex = _Mtx;
 
-    template<class It, class Fn>
-    static IteratorOnly<It, void> for_each(It   from,
-                                           It   to,
-                                           Fn &&fn,
-                                           size_t /* ignore granularity */ = 1)
+    template<class Fn, class It>
+    static IteratorOnly<It, void> loop_(It from, It to, Fn &&fn)
     {
         for (auto it = from; it != to; ++it) fn(*it);
     }
 
-    template<class I, class Fn>
-    static IntegerOnly<I, void> for_each(I    from,
-                                         I    to,
-                                         Fn &&fn,
-                                         size_t /* ignore granularity */ = 1)
+    template<class Fn, class I>
+    static IntegerOnly<I, void> loop_(I from, I to, Fn &&fn)
     {
         for (I i = from; i < to; ++i) fn(i);
     }
+
+    template<class It, class Fn>
+    static void for_each(It   from,
+                         It   to,
+                         Fn &&fn,
+                         size_t /* ignore granularity */ = 1)
+    {
+        loop_(from, to, std::forward<Fn>(fn));
+    }
+
+    template<class I, class MergeFn, class T, class AccessFn>
+    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, access(i)); });
+        return acc;
+    }
+
+    template<class I, class MergeFn, class T>
+    static IteratorOnly<I, T> reduce(I          from,
+                                     I          to,
+                                     const T   &init,
+                                     MergeFn  &&mergefn,
+                                     size_t     /*granularity*/ = 1
+                                     )
+    {
+        return reduce(from, to, init, std::forward<MergeFn>(mergefn),
+                      [](typename I::value_type &i) { return i; });
+    }
 };
 
 using ccr = _ccr<USE_FULL_CONCURRENCY>;
diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp
index 81ef00e6b..937897766 100644
--- a/src/libslic3r/SLA/Rotfinder.cpp
+++ b/src/libslic3r/SLA/Rotfinder.cpp
@@ -1,35 +1,259 @@
 #include <limits>
-#include <exception>
 
-#include <libnest2d/optimizers/nlopt/genetic.hpp>
 #include <libslic3r/SLA/Rotfinder.hpp>
-#include <libslic3r/SLA/SupportTree.hpp>
+#include <libslic3r/SLA/Concurrency.hpp>
+
+#include <libslic3r/Optimize/BruteforceOptimizer.hpp>
+
+#include "libslic3r/SLAPrint.hpp"
+#include "libslic3r/PrintConfig.hpp"
+
+#include <libslic3r/Geometry.hpp>
 #include "Model.hpp"
 
-namespace Slic3r {
-namespace sla {
+#include <thread>
 
-std::array<double, 3> find_best_rotation(const ModelObject& modelobj,
-                                         float accuracy,
-                                         std::function<void(unsigned)> statuscb,
-                                         std::function<bool()> stopcond)
+namespace Slic3r { namespace sla {
+
+inline bool is_on_floor(const SLAPrintObject &mo)
 {
-    using libnest2d::opt::Method;
-    using libnest2d::opt::bound;
-    using libnest2d::opt::Optimizer;
-    using libnest2d::opt::TOptimizer;
-    using libnest2d::opt::StopCriteria;
+    auto opt_elevation = mo.config().support_object_elevation.getFloat();
+    auto opt_padaround = mo.config().pad_around_object.getBool();
 
-    static const unsigned MAX_TRIES = 100000;
+    return opt_elevation < EPSILON || opt_padaround;
+}
+
+// Find transformed mesh ground level without copy and with parallel 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 accessfn = [&mesh, &tr] (size_t vi) {
+        return (tr * mesh.its.vertices[vi].template cast<double>()).z();
+    };
+
+    double zmin = std::numeric_limits<double>::max();
+    size_t granularity = vsize / threads;
+    return ccr_par::reduce(size_t(0), vsize, zmin, minfn, accessfn, granularity);
+}
+
+// Get the vertices of a triangle directly in an array of 3 points
+std::array<Vec3d, 3> get_triangle_vertices(const TriangleMesh &mesh,
+                                           size_t              faceidx)
+{
+    const auto &face = mesh.its.indices[faceidx];
+    return {Vec3d{mesh.its.vertices[face(0)].cast<double>()},
+            Vec3d{mesh.its.vertices[face(1)].cast<double>()},
+            Vec3d{mesh.its.vertices[face(2)].cast<double>()}};
+}
+
+std::array<Vec3d, 3> 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 Facestats {
+    Vec3d  normal;
+    double area;
+
+    explicit Facestats(const std::array<Vec3d, 3> &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 Facestats &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;
+}
+
+template<class AccessFn>
+double sum_score(AccessFn &&accessfn, size_t facecount, size_t Nthreads)
+{
+    double initv     = 0.;
+    auto   mergefn   = std::plus<double>{};
+    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) {
+        Facestats fc{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<Vec3d, 3> tri = get_transformed_triangle(mesh, tr, fi);
+        Facestats fc{tri};
+
+        if (tri[0].z() <= zlvl && tri[1].z() <= zlvl && tri[2].z() <= zlvl)
+            return -fc.area * POINTS_PER_UNIT_AREA;
+
+        return get_score(fc);
+    };
+
+    size_t facecount = mesh.its.indices.size();
+    return sum_score(accessfn, facecount, Nthreads) / facecount;
+}
+
+using XYRotation = std::array<double, 2>;
+
+// 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<size_t N, class Fn, class It, class StopCond>
+std::array<double, N> find_min_score(Fn &&fn, It from, It to, StopCond &&stopfn)
+{
+    std::array<double, N> ret;
+
+    double score = std::numeric_limits<double>::max();
+
+    size_t Nthreads = std::thread::hardware_concurrency();
+    size_t dist = std::distance(from, to);
+    std::vector<double> 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<XYRotation> 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<double>(1e3) * scaled(1.));
+
+    size_t facecount = chull.its.indices.size();
+
+    struct RotArea { XYRotation rot; double area; };
+
+    auto inputs = reserve_vector<RotArea>(facecount);
+
+    auto rotcmp = [](const RotArea &r1, const RotArea &r2) {
+        double xdiff = r1.rot[X] - r2.rot[X], ydiff = r1.rot[Y] - r2.rot[Y];
+        return std::abs(xdiff) < EPSILON ? ydiff < 0. : xdiff < 0.;
+    };
+
+    auto eqcmp = [](const XYRotation &r1, const XYRotation &r2) {
+        double xdiff = r1[X] - r2[X], ydiff = r1[Y] - r2[Y];
+        return std::abs(xdiff) < EPSILON  && std::abs(ydiff) < EPSILON;
+    };
+
+    for (size_t fi = 0; fi < facecount; ++fi) {
+        Facestats fc{get_triangle_vertices(chull, fi)};
+
+        if (fc.area > area_threshold)  {
+            auto q = Eigen::Quaterniond{}.FromTwoVectors(fc.normal, DOWN);
+            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);
+        }
+    }
+
+    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<XYRotation>(std::min(max_count, inputs.size()));
+    for (const RotArea &ra : inputs) ret.emplace_back(ra.rot);
+
+    return ret;
+}
+
+Vec2d find_best_rotation(const SLAPrintObject &        po,
+                         float                         accuracy,
+                         std::function<void(unsigned)> statuscb,
+                         std::function<bool()>         stopcond)
+{
+    static const unsigned MAX_TRIES = 1000;
 
     // return value
-    std::array<double, 3> rot;
+    XYRotation rot;
 
     // We will use only one instance of this converted mesh to examine different
     // rotations
-    const 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
@@ -38,77 +262,61 @@ std::array<double, 3> 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]
-            (double rx, double ry, double rz)
-    {
-        const TriangleMesh& m = mesh;
-
-        // prepare the rotation transformation
-        Transform3d rt = Transform3d::Identity();
-
-        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<double>();
-
-            // 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()));
-        }
-
+    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 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<Method::G_GENETIC> solver(stc);
+    // Different search methods have to be used depending on the model elevation
+    if (is_on_floor(po)) {
 
-    // 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);
+        std::vector<XYRotation> inputs = get_chull_rotations(mesh, max_tries);
+        max_tries = inputs.size();
 
-    // 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);
+        // If the model can be placed on the bed directly, we only need to
+        // check the 3D convex hull face rotations.
 
-    // 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);
+        auto objfn = [&mesh, &statusfn](const XYRotation &rot) {
+            statusfn();
+            Transform3d tr = to_transform3d(rot);
+            return get_model_supportedness_onfloor(mesh, tr);
+        };
 
-    return rot;
+        rot = find_min_score<2>(objfn, inputs.begin(), inputs.end(), stopcond);
+    } else {
+        // Preparing the optimizer.
+        size_t gridsize = std::sqrt(max_tries); // 2D grid has gridsize^2 calls
+        opt::Optimizer<opt::AlgBruteForce> solver(opt::StopCriteria{}
+                                                    .max_iterations(max_tries)
+                                                    .stop_condition(stopcond),
+                                                  gridsize);
+
+        // We are searching rotations around only two axes x, y. Thus the
+        // problem becomes a 2 dimensional optimization task.
+        // We can specify the bounds for a dimension in the following way:
+        auto bounds = opt::bounds({ {-PI, PI}, {-PI, PI} });
+
+        auto result = solver.to_min().optimize(
+            [&mesh, &statusfn] (const XYRotation &rot)
+            {
+                statusfn();
+                return get_model_supportedness(mesh, to_transform3d(rot));
+            }, opt::initvals({0., 0.}), bounds);
+
+        // Save the result and fck off
+        rot = result.optimum;
+    }
+
+    return {rot[0], rot[1]};
 }
 
+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 4469f9731..96561a890 100644
--- a/src/libslic3r/SLA/Rotfinder.hpp
+++ b/src/libslic3r/SLA/Rotfinder.hpp
@@ -4,9 +4,11 @@
 #include <functional>
 #include <array>
 
+#include <libslic3r/Point.hpp>
+
 namespace Slic3r {
 
-class ModelObject;
+class SLAPrintObject;
 
 namespace sla {
 
@@ -25,14 +27,17 @@ namespace sla {
   *
   * @return Returns the rotations around each axis (x, y, z)
   */
-std::array<double, 3> find_best_rotation(
-        const ModelObject& modelobj,
+Vec2d find_best_rotation(
+        const SLAPrintObject& modelobj,
         float accuracy = 1.0f,
         std::function<void(unsigned)> statuscb = [] (unsigned) {},
         std::function<bool()> stopcond = [] () { return false; }
         );
 
-}
-}
+double get_model_supportedness(const SLAPrintObject &mesh,
+                               const Transform3d & tr);
+
+} // namespace sla
+} // namespace Slic3r
 
 #endif // SLAROTFINDER_HPP
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 <libslic3r/SLA/SupportTreeBuildsteps.hpp>
 
 #include <libslic3r/SLA/SpatIndex.hpp>
-#include <libslic3r/Optimizer.hpp>
+#include <libslic3r/Optimize/NLoptOptimizer.hpp>
 #include <boost/log/trivial.hpp>
 
 namespace Slic3r {
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp
index f3e565686..77366c633 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 <GL/glew.h>
 
+#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);
+
+    static constexpr const char * button_txt = L("Optimize orientation");
+    static constexpr const char * slider_txt = L("Accuracy");
+
+    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<std::string> 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;
+
+// TODO:
+
+//    m_rotoptimizewin_state.mobj = ?;
+//    RotoptimzeWindow popup{m_imgui, m_rotoptimizewin_state, {x, y, bottom_limit}};
+
+}
+
 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..c418c4b31 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,79 @@ 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 &          state,
+                         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<std::string>&
     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)) {
diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp
index c847c84b4..978ccf8fc 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"
 
@@ -12,26 +13,41 @@ 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 || int(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)];
 
-    auto r = sla::find_best_rotation(
-        *o,
-        .005f,
+    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;
+//    }
+
+    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
 
     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);
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 <catch2/catch.hpp>
+#include <test_utils.hpp>
+
+#include <libslic3r/Optimize/BruteforceOptimizer.hpp>
+
+#include <libslic3r/Optimize/NLoptOptimizer.hpp>
+
+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<class Opt> 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<class Opt> 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<AlgBruteForce> 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 <libslic3r/SLA/SupportTreeMesher.hpp>
+#include <libslic3r/SLA/Concurrency.hpp>
 
 namespace {
 
@@ -239,3 +240,14 @@ TEST_CASE("halfcone test", "[halfcone]") {
     m.require_shared_vertices();
     m.WriteOBJFile("Halfcone.obj");
 }
+
+TEST_CASE("Test concurrency")
+{
+    std::vector<double> 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<double>{});
+
+    REQUIRE(s == Approx(ref));
+}