diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 09f75c747..263920ecb 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -215,7 +215,8 @@ add_library(libslic3r STATIC
     SimplifyMeshImpl.hpp
     SimplifyMesh.cpp
     MarchingSquares.hpp
-    Optimizer.hpp
+    Optimize/Optimizer.hpp
+    Optimize/NLoptOptimizer.hpp
     ${OpenVDBUtils_SOURCES}
     SLA/Pad.hpp
     SLA/Pad.cpp
diff --git a/src/libslic3r/Optimize/BruteforceOptimizer.hpp b/src/libslic3r/Optimize/BruteforceOptimizer.hpp
new file mode 100644
index 000000000..da4472568
--- /dev/null
+++ b/src/libslic3r/Optimize/BruteforceOptimizer.hpp
@@ -0,0 +1,120 @@
+#ifndef BRUTEFORCEOPTIMIZER_HPP
+#define BRUTEFORCEOPTIMIZER_HPP
+
+#include <libslic3r/Optimize/NLoptOptimizer.hpp>
+
+namespace Slic3r { namespace opt {
+
+namespace detail {
+// Implementing a bruteforce optimizer
+
+template<size_t N>
+constexpr 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;
+}
+
+struct AlgBurteForce {
+    bool to_min;
+    StopCriteria stc;
+    size_t gridsz;
+
+    AlgBurteForce(const StopCriteria &cr, size_t gs): stc{cr}, gridsz{gs} {}
+
+    template<int D, size_t N, class Fn, class Cmp>
+    void run(std::array<size_t, N> &idx,
+             Result<N> &result,
+             const Bounds<N> &bounds,
+             Fn &&fn,
+             Cmp &&cmp)
+    {
+        if (stc.stop_condition()) return;
+
+        if constexpr (D < 0) {
+            Input<N> inp;
+
+            auto max_iter = stc.max_iterations();
+            if (max_iter && num_iter(idx, gridsz) >= max_iter) return;
+
+            for (size_t d = 0; d < N; ++d) {
+                const Bound &b = bounds[d];
+                double step = (b.max() - b.min()) / (gridsz - 1);
+                inp[d] = b.min() + idx[d] * step;
+            }
+
+            auto score = fn(inp);
+            if (cmp(score, result.score)) {
+                result.score = score;
+                result.optimum = inp;
+            }
+
+        } else {
+            for (size_t i = 0; i < gridsz; ++i) {
+                idx[D] = i;
+                run<D - 1>(idx, result, bounds, std::forward<Fn>(fn),
+                           std::forward<Cmp>(cmp));
+            }
+        }
+    }
+
+    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 bruteforce_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 1c94f3c1e..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,24 +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>;
 using AlgNLoptDIRECT  = detail::NLoptAlg<NLOPT_GN_DIRECT>;
-
-// 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 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..f82d6a39e 100644
--- a/src/libslic3r/SLA/Concurrency.hpp
+++ b/src/libslic3r/SLA/Concurrency.hpp
@@ -4,6 +4,7 @@
 #include <tbb/spin_mutex.h>
 #include <tbb/mutex.h>
 #include <tbb/parallel_for.h>
+#include <tbb/parallel_reduce.h>
 #include <algorithm>
 #include <libslic3r/libslic3r.h>
 
@@ -21,28 +22,43 @@ 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 Fn, class MergeFn, class T>
+    static T reduce(I         from,
+                    I         to,
+                    const T & init,
+                    Fn &&     fn,
+                    MergeFn &&mergefn,
+                    size_t    granularity = 1)
     {
-        tbb::parallel_for(tbb::blocked_range{from, to, granularity},
-                          [&fn](const auto &range) {
-            for (I i = range.begin(); i < range.end(); ++i) fn(i);
-        });
+        return tbb::parallel_reduce(
+            tbb::blocked_range{from, to, granularity}, init,
+            [&](const auto &range, T subinit) {
+                T acc = subinit;
+                loop_(range, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); });
+                return acc;
+            },
+            std::forward<MergeFn>(mergefn));
     }
 };
 
@@ -55,23 +71,39 @@ 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 Fn, class MergeFn, class T>
+    static IntegerOnly<I, T> reduce(I         from,
+                                    I         to,
+                                    const T & init,
+                                    Fn &&     fn,
+                                    MergeFn &&mergefn,
+                                    size_t    /*granularity*/ = 1)
+    {
+        T acc = init;
+        loop_(from, to, [&](auto &i) { acc = mergefn(acc, fn(i, acc)); });
+        return acc;
+    }
 };
 
 using ccr = _ccr<USE_FULL_CONCURRENCY>;
diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp
index b4b1fae39..723e50eeb 100644
--- a/src/libslic3r/SLA/Rotfinder.cpp
+++ b/src/libslic3r/SLA/Rotfinder.cpp
@@ -2,23 +2,19 @@
 #include <exception>
 
 //#include <libnest2d/optimizers/nlopt/genetic.hpp>
-#include <libslic3r/Optimizer.hpp>
+#include <libslic3r/Optimize/BruteforceOptimizer.hpp>
 #include <libslic3r/SLA/Rotfinder.hpp>
+#include <libslic3r/SLA/Concurrency.hpp>
 #include <libslic3r/SLA/SupportTree.hpp>
 #include <libslic3r/SLA/SupportPointGenerator.hpp>
 #include <libslic3r/SimplifyMesh.hpp>
 #include "Model.hpp"
 
+#include <thread>
+
 namespace Slic3r {
 namespace sla {
 
-double area(const Vec3d &p1, const Vec3d &p2, const Vec3d &p3) {
-    Vec3d a = p2 - p1;
-    Vec3d b = p3 - p1;
-    Vec3d c = a.cross(b);
-    return 0.5 * c.norm();
-}
-
 using VertexFaceMap = std::vector<std::vector<size_t>>;
 
 VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) {
@@ -35,61 +31,75 @@ VertexFaceMap create_vertex_face_map(const TriangleMesh &mesh) {
     return vmap;
 }
 
+// Find transformed mesh ground level without copy and with parallell reduce.
+double find_ground_level(const TriangleMesh &mesh,
+                         const Transform3d & tr,
+                         size_t              threads)
+{
+    size_t vsize = mesh.its.vertices.size();
+
+    auto minfn = [](double a, double b) { return std::min(a, b); };
+
+    auto findminz = [&mesh, &tr] (size_t vi, double submin) {
+        Vec3d v = tr * mesh.its.vertices[vi].template cast<double>();
+        return std::min(submin, v.z());
+    };
+
+    double zmin = mesh.its.vertices.front().z();
+
+    return ccr_par::reduce(size_t(0), vsize, zmin, findminz, minfn,
+                           vsize / threads);
+}
+
 // Try to guess the number of support points needed to support a mesh
 double calculate_model_supportedness(const TriangleMesh & mesh,
-                                     const VertexFaceMap &vmap,
+//                                     const VertexFaceMap &vmap,
                                      const Transform3d &  tr)
 {
-    static const double POINTS_PER_UNIT_AREA = 1.;
-    static const Vec3d DOWN = {0., 0., -1.};
+    static constexpr double POINTS_PER_UNIT_AREA = 1.;
 
-    double score = 0.;
+    if (mesh.its.vertices.empty()) return std::nan("");
 
-//    double zmin = mesh.bounding_box().min.z();
+    size_t Nthr = std::thread::hardware_concurrency();
+    size_t facesize = mesh.its.indices.size();
 
-//    std::vector<Vec3d> normals(mesh.its.indices.size(), Vec3d::Zero());
+    double zmin = find_ground_level(mesh, tr, Nthr);
 
-    double zmin = 0;
-    for (auto & v : mesh.its.vertices)
-        zmin = std::min(zmin, double((tr * v.cast<double>()).z()));
+    auto score_mergefn = [&mesh, &tr, zmin](size_t fi, double subscore) {
+
+        static const Vec3d DOWN = {0., 0., -1.};
 
-    for (size_t fi = 0; fi < mesh.its.indices.size(); ++fi) {
         const auto &face = mesh.its.indices[fi];
-        Vec3d p1 = tr * mesh.its.vertices[face(0)].cast<double>();
-        Vec3d p2 = tr * mesh.its.vertices[face(1)].cast<double>();
-        Vec3d p3 = tr * mesh.its.vertices[face(2)].cast<double>();
+        Vec3d p1 = tr * mesh.its.vertices[face(0)].template cast<double>();
+        Vec3d p2 = tr * mesh.its.vertices[face(1)].template cast<double>();
+        Vec3d p3 = tr * mesh.its.vertices[face(2)].template cast<double>();
 
-//        auto triang = std::array<Vec3d, 3> {p1, p2, p3};
-//        double a = area(triang.begin(), triang.end());
-        double a = area(p1, p2, p3);
+        Vec3d U = p2 - p1;
+        Vec3d V = p3 - p1;
+        Vec3d C = U.cross(V);
+        Vec3d N = C.normalized();
+        double area = 0.5 * C.norm();
 
         double zlvl = zmin + 0.1;
         if (p1.z() <= zlvl && p2.z() <= zlvl && p3.z() <= zlvl) {
-            score += a * POINTS_PER_UNIT_AREA;
-            continue;
+            //                score += area * POINTS_PER_UNIT_AREA;
+            return subscore;
         }
 
+        double phi = 1. - std::acos(N.dot(DOWN)) / PI;
+        phi = phi * (phi > 0.5);
 
-        Eigen::Vector3d U = p2 - p1;
-        Eigen::Vector3d V = p3 - p1;
-        Vec3d           N = U.cross(V).normalized();
+        //                    std::cout << "area: " << area << std::endl;
 
-        double phi = std::acos(N.dot(DOWN)) / PI;
+        subscore += area * POINTS_PER_UNIT_AREA * phi;
 
-        std::cout << "area: " << a << std::endl;
+        return subscore;
+    };
 
-        score += a * POINTS_PER_UNIT_AREA * phi;
-//        normals[fi] = N;
-    }
+    double score = ccr_seq::reduce(size_t(0), facesize, 0., score_mergefn,
+                                   std::plus<double>{}, facesize / Nthr);
 
-//    for (size_t vi = 0; vi < mesh.its.vertices.size(); ++vi) {
-//        const std::vector<size_t> &neighbors = vmap[vi];
-
-//        const auto &v = mesh.its.vertices[vi];
-//        Vec3d vt =  tr * v.cast<double>();
-//    }
-
-    return score;
+    return score / mesh.its.indices.size();
 }
 
 std::array<double, 2> find_best_rotation(const ModelObject& modelobj,
@@ -97,7 +107,7 @@ std::array<double, 2> find_best_rotation(const ModelObject& modelobj,
                                          std::function<void(unsigned)> statuscb,
                                          std::function<bool()> stopcond)
 {
-    static const unsigned MAX_TRIES = 1000000;
+    static const unsigned MAX_TRIES = 100;
 
     // return value
     std::array<double, 2> rot;
@@ -126,12 +136,14 @@ std::array<double, 2> find_best_rotation(const ModelObject& modelobj,
     auto objfunc = [&mesh, &status, &statuscb, &stopcond, max_tries]
         (const opt::Input<2> &in)
     {
+        std::cout << "in: " << in[0] << " " << in[1] << std::endl;
+
         // prepare the rotation transformation
         Transform3d rt = Transform3d::Identity();
         rt.rotate(Eigen::AngleAxisd(in[1], Vec3d::UnitY()));
         rt.rotate(Eigen::AngleAxisd(in[0], Vec3d::UnitX()));
 
-        double score = sla::calculate_model_supportedness(mesh, {}, rt);
+        double score = sla::calculate_model_supportedness(mesh, rt);
         std::cout << score << std::endl;
 
         // report status
@@ -142,10 +154,11 @@ std::array<double, 2> find_best_rotation(const ModelObject& modelobj,
 
     // Firing up the genetic optimizer. For now it uses the nlopt library.
 
-    opt::Optimizer<opt::AlgNLoptDIRECT> solver(opt::StopCriteria{}
-                                           .max_iterations(max_tries)
-                                           .rel_score_diff(1e-3)
-                                           .stop_condition(stopcond));
+    opt::Optimizer<opt::AlgBruteForce> solver(opt::StopCriteria{}
+                                                  .max_iterations(max_tries)
+                                                  .rel_score_diff(1e-6)
+                                                  .stop_condition(stopcond),
+                                              10 /*grid size*/);
 
     // We are searching rotations around the three axes x, y, z. Thus the
     // problem becomes a 3 dimensional optimization task.
@@ -153,7 +166,7 @@ std::array<double, 2> find_best_rotation(const ModelObject& modelobj,
     auto b = opt::Bound{-PI, PI};
 
     // Now we start the optimization process with initial angles (0, 0, 0)
-    auto result = solver.to_max().optimize(objfunc, opt::initvals({0.0, 0.0}),
+    auto result = solver.to_min().optimize(objfunc, opt::initvals({0.0, 0.0}),
                                            opt::bounds({b, b}));
 
     // Save the result and fck off
diff --git a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
index 0e7af8d50..3c39c64e6 100644
--- a/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
+++ b/src/libslic3r/SLA/SupportTreeBuildsteps.cpp
@@ -1,7 +1,7 @@
 #include <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 {