diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 4a762f7e1..2abe94656 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -219,6 +219,9 @@ add_library(libslic3r STATIC
     SimplifyMeshImpl.hpp
     SimplifyMesh.cpp
     MarchingSquares.hpp
+    Execution/Execution.hpp
+    Execution/ExecutionSeq.hpp
+    Execution/ExecutionTBB.hpp
     Optimize/Optimizer.hpp
     Optimize/NLoptOptimizer.hpp
     Optimize/BruteforceOptimizer.hpp
diff --git a/src/libslic3r/Execution/Execution.hpp b/src/libslic3r/Execution/Execution.hpp
new file mode 100644
index 000000000..e4bad9f23
--- /dev/null
+++ b/src/libslic3r/Execution/Execution.hpp
@@ -0,0 +1,132 @@
+#ifndef EXECUTION_HPP
+#define EXECUTION_HPP
+
+#include <type_traits>
+#include <utility>
+#include <cstddef>
+#include <iterator>
+
+#include "libslic3r/libslic3r.h"
+
+namespace Slic3r {
+
+// Borrowed from C++20
+template<class T>
+using remove_cvref_t = std::remove_reference_t<std::remove_cv_t<T>>;
+
+// Override for valid execution policies
+template<class EP> struct IsExecutionPolicy_ : public std::false_type {};
+
+template<class EP> constexpr bool IsExecutionPolicy =
+    IsExecutionPolicy_<remove_cvref_t<EP>>::value;
+
+template<class EP, class T = void>
+using ExecutionPolicyOnly = std::enable_if_t<IsExecutionPolicy<EP>, T>;
+
+namespace execution {
+
+// This struct needs to be specialized for each execution policy.
+// See ExecutionSeq.hpp and ExecutionTBB.hpp for example.
+template<class EP, class En = void> struct Traits {};
+
+template<class EP> using AsTraits = Traits<remove_cvref_t<EP>>;
+
+// Each execution policy should declare two types of mutexes. A a spin lock and
+// a blocking mutex. These types should satisfy the BasicLockable concept.
+template<class EP> using SpinningMutex = typename Traits<EP>::SpinningMutex;
+template<class EP> using BlockingMutex = typename Traits<EP>::BlockingMutex;
+
+// Query the available threads for concurrency.
+template<class EP, class = ExecutionPolicyOnly<EP> >
+size_t max_concurrency(const EP &ep)
+{
+    return AsTraits<EP>::max_concurrency(ep);
+}
+
+// foreach loop with the execution policy passed as argument. Granularity can
+// be specified explicitly. max_concurrency() can be used for optimal results.
+template<class EP, class It, class Fn, class = ExecutionPolicyOnly<EP>>
+void for_each(const EP &ep, It from, It to, Fn &&fn, size_t granularity = 1)
+{
+    AsTraits<EP>::for_each(ep, from, to, std::forward<Fn>(fn), granularity);
+}
+
+// A reduce operation with the execution policy passed as argument.
+// mergefn has T(const T&, const T&) signature
+// accessfn has T(I) signature if I is an integral type and
+// T(const I::value_type &) if I is an iterator type.
+template<class EP,
+         class I,
+         class MergeFn,
+         class T,
+         class AccessFn,
+         class = ExecutionPolicyOnly<EP> >
+T reduce(const EP & ep,
+         I          from,
+         I          to,
+         const T &  init,
+         MergeFn && mergefn,
+         AccessFn &&accessfn,
+         size_t     granularity = 1)
+{
+    return AsTraits<EP>::reduce(ep, from, to, init,
+                                std::forward<MergeFn>(mergefn),
+                                std::forward<AccessFn>(accessfn),
+                                granularity);
+}
+
+// An overload of reduce method to be used with iterators as 'from' and 'to'
+// arguments. Access functor is omitted here.
+template<class EP,
+         class I,
+         class MergeFn,
+         class T,
+         class = ExecutionPolicyOnly<EP> >
+T reduce(const EP &ep,
+         I         from,
+         I         to,
+         const T & init,
+         MergeFn &&mergefn,
+         size_t    granularity = 1)
+{
+    return reduce(
+        ep, from, to, init, std::forward<MergeFn>(mergefn),
+        [](const auto &i) { return i; }, granularity);
+}
+
+template<class EP,
+         class I,
+         class T,
+         class AccessFn,
+         class = ExecutionPolicyOnly<EP>>
+T accumulate(const EP & ep,
+             I          from,
+             I          to,
+             const T &  init,
+             AccessFn &&accessfn,
+             size_t     granularity = 1)
+{
+    return reduce(ep, from, to, init, std::plus<T>{},
+                  std::forward<AccessFn>(accessfn), granularity);
+}
+
+
+template<class EP,
+         class I,
+         class T,
+         class = ExecutionPolicyOnly<EP> >
+T accumulate(const EP &ep,
+             I         from,
+             I         to,
+             const T & init,
+             size_t    granularity = 1)
+{
+    return reduce(
+        ep, from, to, init, std::plus<T>{}, [](const auto &i) { return i; },
+        granularity);
+}
+
+} // namespace execution_policy
+} // namespace Slic3r
+
+#endif // EXECUTION_HPP
diff --git a/src/libslic3r/Execution/ExecutionSeq.hpp b/src/libslic3r/Execution/ExecutionSeq.hpp
new file mode 100644
index 000000000..321d65631
--- /dev/null
+++ b/src/libslic3r/Execution/ExecutionSeq.hpp
@@ -0,0 +1,84 @@
+#ifndef EXECUTIONSEQ_HPP
+#define EXECUTIONSEQ_HPP
+
+#ifdef PRUSASLICER_USE_EXECUTION_STD // Conflicts with our version of TBB
+#include <execution>
+#endif
+
+#include "Execution.hpp"
+
+namespace Slic3r {
+
+// Execution policy implementing dummy sequential algorithms
+struct ExecutionSeq {};
+
+template<> struct IsExecutionPolicy_<ExecutionSeq> : public std::true_type {};
+
+static constexpr ExecutionSeq ex_seq = {};
+
+template<class EP> struct IsSequentialEP_ { static constexpr bool value = false; };
+
+template<> struct IsSequentialEP_<ExecutionSeq>: public std::true_type {};
+#ifdef PRUSASLICER_USE_EXECUTION_STD
+template<> struct IsExecutionPolicy_<std::execution::sequenced_policy>: public std::true_type {};
+template<> struct IsSequentialEP_<std::execution::sequenced_policy>: public std::true_type {};
+#endif
+
+template<class EP>
+constexpr bool IsSequentialEP = IsSequentialEP_<remove_cvref_t<EP>>::value;
+
+template<class EP, class R = EP>
+using SequentialEPOnly = std::enable_if_t<IsSequentialEP<EP>, R>;
+
+template<class EP>
+struct execution::Traits<EP, SequentialEPOnly<EP, void>> {
+private:
+    struct _Mtx { inline void lock() {} inline void unlock() {} };
+
+    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 Fn, class I>
+    static IntegerOnly<I, void> loop_(I from, I to, Fn &&fn)
+    {
+        for (I i = from; i < to; ++i) fn(i);
+    }
+
+public:
+    using SpinningMutex = _Mtx;
+    using BlockingMutex = _Mtx;
+
+    template<class It, class Fn>
+    static void for_each(const EP &,
+                         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(const EP &,
+                    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;
+    }
+
+    static size_t max_concurrency(const EP &) { return 1; }
+};
+
+} // namespace Slic3r
+
+#endif // EXECUTIONSEQ_HPP
diff --git a/src/libslic3r/Execution/ExecutionTBB.hpp b/src/libslic3r/Execution/ExecutionTBB.hpp
new file mode 100644
index 000000000..cf6373c46
--- /dev/null
+++ b/src/libslic3r/Execution/ExecutionTBB.hpp
@@ -0,0 +1,77 @@
+#ifndef EXECUTIONTBB_HPP
+#define EXECUTIONTBB_HPP
+
+#include <tbb/spin_mutex.h>
+#include <tbb/mutex.h>
+#include <tbb/parallel_for.h>
+#include <tbb/parallel_reduce.h>
+#include <tbb/task_arena.h>
+
+#include "Execution.hpp"
+
+namespace Slic3r {
+
+struct ExecutionTBB {};
+template<> struct IsExecutionPolicy_<ExecutionTBB> : public std::true_type {};
+
+// Execution policy using Intel TBB library under the hood.
+static constexpr ExecutionTBB ex_tbb = {};
+
+template<> struct execution::Traits<ExecutionTBB> {
+private:
+
+    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);
+    }
+
+public:
+    using SpinningMutex = tbb::spin_mutex;
+    using BlockingMutex = tbb::mutex;
+
+    template<class It, class Fn>
+    static void for_each(const ExecutionTBB &,
+                         It from, It to, Fn &&fn, size_t granularity)
+    {
+        tbb::parallel_for(tbb::blocked_range{from, to, granularity},
+                          [&fn](const auto &range) {
+            loop_(range, std::forward<Fn>(fn));
+        });
+    }
+
+    template<class I, class MergeFn, class T, class AccessFn>
+    static T reduce(const ExecutionTBB &,
+                    I          from,
+                    I          to,
+                    const T   &init,
+                    MergeFn  &&mergefn,
+                    AccessFn &&access,
+                    size_t     granularity = 1
+                    )
+    {
+        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));
+    }
+
+    static size_t max_concurrency(const ExecutionTBB &)
+    {
+        return tbb::this_task_arena::max_concurrency();
+    }
+};
+
+}
+
+#endif // EXECUTIONTBB_HPP
diff --git a/src/libslic3r/MTUtils.hpp b/src/libslic3r/MTUtils.hpp
index 555cfe501..7b903f66c 100644
--- a/src/libslic3r/MTUtils.hpp
+++ b/src/libslic3r/MTUtils.hpp
@@ -106,13 +106,8 @@ template<class C> bool all_of(const C &container)
                        });
 }
 
-template<class T> struct remove_cvref
-{
-    using type =
-        typename std::remove_cv<typename std::remove_reference<T>::type>::type;
-};
-
-template<class T> using remove_cvref_t = typename remove_cvref<T>::type;
+template<class T>
+using remove_cvref_t = std::remove_reference_t<std::remove_cv_t<T>>;
 
 /// Exactly like Matlab https://www.mathworks.com/help/matlab/ref/linspace.html
 template<class T, class I, class = IntegerOnly<I>>
diff --git a/src/libslic3r/SLA/Concurrency.hpp b/src/libslic3r/SLA/Concurrency.hpp
index 8ff0ff809..7299101b3 100644
--- a/src/libslic3r/SLA/Concurrency.hpp
+++ b/src/libslic3r/SLA/Concurrency.hpp
@@ -1,16 +1,10 @@
 #ifndef SLA_CONCURRENCY_H
 #define SLA_CONCURRENCY_H
 
-#include <tbb/spin_mutex.h>
-#include <tbb/mutex.h>
-#include <tbb/parallel_for.h>
-#include <tbb/parallel_reduce.h>
-#include <tbb/task_arena.h>
+// FIXME: Deprecated
 
-#include <algorithm>
-#include <numeric>
-
-#include <libslic3r/libslic3r.h>
+#include <libslic3r/Execution/ExecutionSeq.hpp>
+#include <libslic3r/Execution/ExecutionTBB.hpp>
 
 namespace Slic3r {
 namespace sla {
@@ -23,124 +17,48 @@ template<bool> struct _ccr {};
 
 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);
-    }
+    using SpinningMutex = execution::SpinningMutex<ExecutionTBB>;
+    using BlockingMutex = execution::BlockingMutex<ExecutionTBB>;
 
     template<class It, class Fn>
     static void for_each(It from, It to, Fn &&fn, size_t granularity = 1)
     {
-        tbb::parallel_for(tbb::blocked_range{from, to, granularity},
-                          [&fn](const auto &range) {
-            loop_(range, std::forward<Fn>(fn));
-        });
+        execution::for_each(ex_tbb, from, to, std::forward<Fn>(fn), granularity);
     }
 
-    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
-                    )
+    template<class...Args>
+    static auto reduce(Args&&...args)
     {
-        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);
+        return execution::reduce(ex_tbb, std::forward<Args>(args)...);
     }
 
     static size_t max_concurreny()
     {
-        return tbb::this_task_arena::max_concurrency();
+        return execution::max_concurrency(ex_tbb);
     }
 };
 
 template<> struct _ccr<false>
 {
-private:
-    struct _Mtx { inline void lock() {} inline void unlock() {} };
-
-public:
-    using SpinningMutex = _Mtx;
-    using BlockingMutex = _Mtx;
-
-    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 Fn, class I>
-    static IntegerOnly<I, void> loop_(I from, I to, Fn &&fn)
-    {
-        for (I i = from; i < to; ++i) fn(i);
-    }
+    using SpinningMutex = execution::SpinningMutex<ExecutionSeq>;
+    using BlockingMutex = execution::BlockingMutex<ExecutionSeq>;
 
     template<class It, class Fn>
-    static void for_each(It   from,
-                         It   to,
-                         Fn &&fn,
-                         size_t /* ignore granularity */ = 1)
+    static void for_each(It from, It to, Fn &&fn, size_t granularity = 1)
     {
-        loop_(from, to, std::forward<Fn>(fn));
+        execution::for_each(ex_seq, from, to, std::forward<Fn>(fn), granularity);
     }
 
-    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
-                    )
+    template<class...Args>
+    static auto reduce(Args&&...args)
     {
-        T acc = init;
-        loop_(from, to, [&](auto &i) { acc = mergefn(acc, access(i)); });
-        return acc;
+        return execution::reduce(ex_seq, std::forward<Args>(args)...);
     }
 
-    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
-                                     )
+    static size_t max_concurreny()
     {
-        return reduce(from, to, init, std::forward<MergeFn>(mergefn),
-                      [](typename I::value_type &i) { return i; });
+        return execution::max_concurrency(ex_seq);
     }
-
-    static size_t max_concurreny() { return 1; }
 };
 
 using ccr = _ccr<USE_FULL_CONCURRENCY>;
diff --git a/src/libslic3r/SLA/Rotfinder.cpp b/src/libslic3r/SLA/Rotfinder.cpp
index 702690c19..b84921279 100644
--- a/src/libslic3r/SLA/Rotfinder.cpp
+++ b/src/libslic3r/SLA/Rotfinder.cpp
@@ -1,91 +1,110 @@
 #include <limits>
 
 #include <libslic3r/SLA/Rotfinder.hpp>
-#include <libslic3r/SLA/Concurrency.hpp>
+
+#include <libslic3r/Execution/ExecutionTBB.hpp>
+#include <libslic3r/Execution/ExecutionSeq.hpp>
 
 #include <libslic3r/Optimize/BruteforceOptimizer.hpp>
+#include <libslic3r/Optimize/NLoptOptimizer.hpp>
 
 #include "libslic3r/SLAPrint.hpp"
 #include "libslic3r/PrintConfig.hpp"
 
 #include <libslic3r/Geometry.hpp>
-#include "Model.hpp"
 
 #include <thread>
 
 namespace Slic3r { namespace sla {
 
-inline bool is_on_floor(const SLAPrintObject &mo)
-{
-    auto opt_elevation = mo.config().support_object_elevation.getFloat();
-    auto opt_padaround = mo.config().pad_around_object.getBool();
+namespace {
 
-    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);
-}
+inline const Vec3f DOWN = {0.f, 0.f, -1.f};
+constexpr double POINTS_PER_UNIT_AREA = 1.f;
 
 // Get the vertices of a triangle directly in an array of 3 points
-std::array<Vec3d, 3> get_triangle_vertices(const TriangleMesh &mesh,
+std::array<Vec3f, 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>()}};
+    return {mesh.its.vertices[face(0)],
+            mesh.its.vertices[face(1)],
+            mesh.its.vertices[face(2)]};
 }
 
-std::array<Vec3d, 3> get_transformed_triangle(const TriangleMesh &mesh,
-                                              const Transform3d & tr,
+std::array<Vec3f, 3> get_transformed_triangle(const TriangleMesh &mesh,
+                                              const Transform3f & tr,
                                               size_t              faceidx)
 {
     const auto &tri = get_triangle_vertices(mesh, faceidx);
     return {tr * tri[0], tr * tri[1], tr * tri[2]};
 }
 
+template<class T> Vec<3, T> normal(const std::array<Vec<3, T>, 3> &tri)
+{
+    Vec<3, T> U = tri[1] - tri[0];
+    Vec<3, T> V = tri[2] - tri[0];
+    return U.cross(V).normalized();
+}
+
+template<class T, class AccessFn>
+T sum_score(AccessFn &&accessfn, size_t facecount, size_t Nthreads)
+{
+    T  initv         = 0.;
+    auto   mergefn   = [](T a, T b) { return a + b; };
+    size_t grainsize = facecount / Nthreads;
+    size_t from = 0, to = facecount;
+
+    return execution::reduce(ex_tbb, from, to, initv, mergefn, accessfn, grainsize);
+}
+
+// Try to guess the number of support points needed to support a mesh
+double get_misalginment_score(const TriangleMesh &mesh, const Transform3f &tr)
+{
+    if (mesh.its.vertices.empty()) return std::nan("");
+
+    auto accessfn = [&mesh, &tr](size_t fi) {
+        auto triangle = get_transformed_triangle(mesh, tr, fi);
+        Vec3f U = triangle[1] - triangle[0];
+        Vec3f V = triangle[2] - triangle[0];
+        Vec3f C = U.cross(V);
+
+        // We should score against the alignment with the reference planes
+        return scaled<int_fast64_t>(std::abs(C.dot(Vec3f::UnitX())) +
+                                    std::abs(C.dot(Vec3f::UnitY())));
+    };
+
+    size_t facecount = mesh.its.indices.size();
+    size_t Nthreads  = std::thread::hardware_concurrency();
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
+}
+
 // Get area and normal of a triangle
 struct Facestats {
-    Vec3d  normal;
+    Vec3f  normal;
     double area;
 
-    explicit Facestats(const std::array<Vec3d, 3> &triangle)
+    explicit Facestats(const std::array<Vec3f, 3> &triangle)
     {
-        Vec3d U = triangle[1] - triangle[0];
-        Vec3d V = triangle[2] - triangle[0];
-        Vec3d C = U.cross(V);
+        Vec3f U = triangle[1] - triangle[0];
+        Vec3f V = triangle[2] - triangle[0];
+        Vec3f 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)
+inline double get_supportedness_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;
+    float phi = 1. - std::acos(fc.normal.dot(DOWN)) / float(PI);
 
-    // Only consider faces that have have slopes below 90 deg:
-    phi = phi * (phi > 0.5);
+    // Only consider faces that have slopes below 90 deg:
+    phi = phi * (phi >= 0.5f);
 
     // Make the huge slopes more significant than the smaller slopes
     phi = phi * phi * phi;
@@ -94,96 +113,92 @@ inline double get_score(const Facestats &fc)
     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)
+double get_supportedness_score(const TriangleMesh &mesh, const Transform3f &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);
+
+        return get_supportedness_score(fc);
     };
 
     size_t facecount = mesh.its.indices.size();
     size_t Nthreads  = std::thread::hardware_concurrency();
-    return sum_score(accessfn, facecount, Nthreads) / facecount;
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
 }
 
-double get_model_supportedness_onfloor(const TriangleMesh &mesh,
-                                       const Transform3d & tr)
+// Find transformed mesh ground level without copy and with parallel reduce.
+float find_ground_level(const TriangleMesh &mesh,
+                         const Transform3f & tr,
+                         size_t              threads)
+{
+    size_t vsize = mesh.its.vertices.size();
+
+    auto minfn = [](float a, float b) { return std::min(a, b); };
+
+    auto accessfn = [&mesh, &tr] (size_t vi) {
+        return (tr * mesh.its.vertices[vi]).z();
+    };
+
+    auto zmin = std::numeric_limits<float>::max();
+    size_t granularity = vsize / threads;
+    return execution::reduce(ex_tbb, size_t(0), vsize, zmin, minfn, accessfn, granularity);
+}
+
+float get_supportedness_onfloor_score(const TriangleMesh &mesh,
+                                      const Transform3f & 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
+    float zmin = find_ground_level(mesh, tr, Nthreads);
+    float zlvl = zmin + 0.1f; // 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);
+        std::array<Vec3f, 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);
+        return get_supportedness_score(fc);
     };
 
     size_t facecount = mesh.its.indices.size();
-    return sum_score(accessfn, facecount, Nthreads) / facecount;
+    double S = unscaled(sum_score<int_fast64_t>(accessfn, facecount, Nthreads));
+
+    return S / facecount;
 }
 
 using XYRotation = std::array<double, 2>;
 
 // prepare the rotation transformation
-Transform3d to_transform3d(const XYRotation &rot)
+Transform3f to_transform3f(const XYRotation &rot)
 {
-    Transform3d rt = Transform3d::Identity();
-    rt.rotate(Eigen::AngleAxisd(rot[1], Vec3d::UnitY()));
-    rt.rotate(Eigen::AngleAxisd(rot[0], Vec3d::UnitX()));
+    Transform3f rt = Transform3f::Identity();
+    rt.rotate(Eigen::AngleAxisf(float(rot[1]), Vec3f::UnitY()));
+    rt.rotate(Eigen::AngleAxisf(float(rot[0]), Vec3f::UnitX()));
+
     return rt;
 }
 
-XYRotation from_transform3d(const Transform3d &tr)
+XYRotation from_transform3f(const Transform3f &tr)
 {
-    Vec3d rot3d = Geometry::Transformation {tr}.get_rotation();
-    return {rot3d.x(), rot3d.y()};
+    Vec3d rot3 = Geometry::Transformation{tr.cast<double>()}.get_rotation();
+    return {rot3.x(), rot3.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)
+inline bool is_on_floor(const SLAPrintObjectConfig &cfg)
 {
-    std::array<double, N> ret = {};
+    auto opt_elevation = cfg.support_object_elevation.getFloat();
+    auto opt_padaround = cfg.pad_around_object.getBool();
 
-    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;
+    return opt_elevation < EPSILON || opt_padaround;
 }
 
 // collect the rotations for each face of the convex hull
@@ -214,8 +229,8 @@ std::vector<XYRotation> get_chull_rotations(const TriangleMesh &mesh, size_t max
         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);
+            auto q = Eigen::Quaternionf{}.FromTwoVectors(fc.normal, DOWN);
+            XYRotation rot = from_transform3f(Transform3f::Identity() * q);
             RotArea ra = {rot, fc.area};
 
             auto it = std::lower_bound(inputs.begin(), inputs.end(), ra, rotcmp);
@@ -238,10 +253,95 @@ std::vector<XYRotation> get_chull_rotations(const TriangleMesh &mesh, size_t max
     return ret;
 }
 
-Vec2d find_best_rotation(const SLAPrintObject &        po,
-                         float                         accuracy,
-                         std::function<void(unsigned)> statuscb,
-                         std::function<bool()>         stopcond)
+// 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);
+
+    execution::for_each(
+        ex_tbb, 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;
+}
+
+} // namespace
+
+Vec2d find_best_misalignment_rotation(const ModelObject &      mo,
+                                      const RotOptimizeParams &params)
+{
+    static constexpr unsigned MAX_TRIES = 1000;
+
+    // return value
+    XYRotation rot;
+
+    // We will use only one instance of this converted mesh to examine different
+    // rotations
+    TriangleMesh mesh = mo.raw_mesh();
+    mesh.require_shared_vertices();
+
+    // To keep track of the number of iterations
+    int status = 0;
+
+    // The maximum number of iterations
+    auto max_tries = unsigned(params.accuracy() * MAX_TRIES);
+
+    auto &statuscb = params.statuscb();
+
+    // call status callback with zero, because we are at the start
+    statuscb(status);
+
+    auto statusfn = [&statuscb, &status, &max_tries] {
+        // report status
+        statuscb(++status * 100.0/max_tries);
+    };
+
+    auto stopcond = [&statuscb] {
+        return ! statuscb(-1);
+    };
+
+    // Preparing the optimizer.
+    size_t gridsize = std::sqrt(max_tries);
+    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/2, PI/2}, {-PI/2, PI/2} });
+
+    auto result = solver.to_max().optimize(
+        [&mesh, &statusfn] (const XYRotation &rot)
+        {
+            statusfn();
+            return get_misalginment_score(mesh, to_transform3f(rot));
+        }, opt::initvals({0., 0.}), bounds);
+
+    rot = result.optimum;
+
+    return {rot[0], rot[1]};
+}
+
+Vec2d find_least_supports_rotation(const ModelObject &      mo,
+                                   const RotOptimizeParams &params)
 {
     static const unsigned MAX_TRIES = 1000;
 
@@ -250,14 +350,16 @@ Vec2d find_best_rotation(const SLAPrintObject &        po,
 
     // We will use only one instance of this converted mesh to examine different
     // rotations
-    TriangleMesh mesh = po.model_object()->raw_mesh();
+    TriangleMesh mesh = mo.raw_mesh();
     mesh.require_shared_vertices();
 
     // To keep track of the number of iterations
     unsigned status = 0;
 
     // The maximum number of iterations
-    auto max_tries = unsigned(accuracy * MAX_TRIES);
+    auto max_tries = unsigned(params.accuracy() * MAX_TRIES);
+
+    auto &statuscb = params.statuscb();
 
     // call status callback with zero, because we are at the start
     statuscb(status);
@@ -267,8 +369,18 @@ Vec2d find_best_rotation(const SLAPrintObject &        po,
         statuscb(unsigned(++status * 100.0/max_tries) );
     };
 
+    auto stopcond = [&statuscb] {
+        return ! statuscb(-1);
+    };
+
+    SLAPrintObjectConfig pocfg;
+    if (params.print_config())
+        pocfg.apply(*params.print_config(), true);
+
+    pocfg.apply(mo.config.get());
+
     // Different search methods have to be used depending on the model elevation
-    if (is_on_floor(po)) {
+    if (is_on_floor(pocfg)) {
 
         std::vector<XYRotation> inputs = get_chull_rotations(mesh, max_tries);
         max_tries = inputs.size();
@@ -278,17 +390,19 @@ Vec2d find_best_rotation(const SLAPrintObject &        po,
 
         auto objfn = [&mesh, &statusfn](const XYRotation &rot) {
             statusfn();
-            Transform3d tr = to_transform3d(rot);
-            return get_model_supportedness_onfloor(mesh, tr);
+            Transform3f tr = to_transform3f(rot);
+            return get_supportedness_onfloor_score(mesh, tr);
         };
 
         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),
+                                                      .max_iterations(max_tries)
+                                                      .stop_condition(stopcond),
                                                   gridsize);
 
         // We are searching rotations around only two axes x, y. Thus the
@@ -300,23 +414,14 @@ Vec2d find_best_rotation(const SLAPrintObject &        po,
             [&mesh, &statusfn] (const XYRotation &rot)
             {
                 statusfn();
-                return get_model_supportedness(mesh, to_transform3d(rot));
+                return get_supportedness_score(mesh, to_transform3f(rot));
             }, opt::initvals({0., 0.}), bounds);
 
-        // Save the result and fck off
+        // Save the result
         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 96561a890..77a39016d 100644
--- a/src/libslic3r/SLA/Rotfinder.hpp
+++ b/src/libslic3r/SLA/Rotfinder.hpp
@@ -8,10 +8,39 @@
 
 namespace Slic3r {
 
+class ModelObject;
 class SLAPrintObject;
+class TriangleMesh;
+class DynamicPrintConfig;
 
 namespace sla {
 
+using RotOptimizeStatusCB = std::function<bool(int)>;
+
+class RotOptimizeParams {
+    float m_accuracy = 1.;
+    const DynamicPrintConfig *m_print_config = nullptr;
+    RotOptimizeStatusCB m_statuscb = [](int) { return true; };
+
+public:
+
+    RotOptimizeParams &accuracy(float a) { m_accuracy = a; return *this; }
+    RotOptimizeParams &print_config(const DynamicPrintConfig *c)
+    {
+        m_print_config = c;
+        return *this;
+    }
+    RotOptimizeParams &statucb(RotOptimizeStatusCB cb)
+    {
+        m_statuscb = std::move(cb);
+        return *this;
+    }
+
+    float accuracy() const { return m_accuracy; }
+    const DynamicPrintConfig * print_config() const { return m_print_config; }
+    const RotOptimizeStatusCB &statuscb() const { return m_statuscb; }
+};
+
 /**
   * The function should find the best rotation for SLA upside down printing.
   *
@@ -19,23 +48,22 @@ namespace sla {
   * @param accuracy The optimization accuracy from 0.0f to 1.0f. Currently,
   * the nlopt genetic optimizer is used and the number of iterations is
   * accuracy * 100000. This can change in the future.
-  * @param statuscb A status indicator callback called with the unsigned
+  * @param statuscb A status indicator callback called with the int
   * argument spanning from 0 to 100. May not reach 100 if the optimization finds
-  * an optimum before max iterations are reached.
-  * @param stopcond A function that if returns true, the search process will be
-  * terminated and the best solution found will be returned.
+  * an optimum before max iterations are reached. It should return a boolean
+  * signaling if the operation may continue (true) or not (false). A status
+  * value lower than 0 shall not update the status but still return a valid
+  * continuation indicator.
   *
   * @return Returns the rotations around each axis (x, y, z)
   */
-Vec2d find_best_rotation(
-        const SLAPrintObject& modelobj,
-        float accuracy = 1.0f,
-        std::function<void(unsigned)> statuscb = [] (unsigned) {},
-        std::function<bool()> stopcond = [] () { return false; }
-        );
+Vec2d find_best_misalignment_rotation(const ModelObject &modelobj,
+                                      const RotOptimizeParams & = {});
 
-double get_model_supportedness(const SLAPrintObject &mesh,
-                               const Transform3d & tr);
+Vec2d find_least_supports_rotation(const ModelObject &modelobj,
+                                   const RotOptimizeParams & = {});
+
+double find_Z_fit_to_bed_rotation(const ModelObject &mo, const BoundingBox &bed);
 
 } // namespace sla
 } // namespace Slic3r
diff --git a/src/slic3r/GUI/GUI_Factories.cpp b/src/slic3r/GUI/GUI_Factories.cpp
index 0bd523b26..01ea53ad3 100644
--- a/src/slic3r/GUI/GUI_Factories.cpp
+++ b/src/slic3r/GUI/GUI_Factories.cpp
@@ -858,10 +858,6 @@ void MenuFactory::create_sla_object_menu()
         []() { return plater()->can_split(true); }, m_parent);
 
     m_sla_object_menu.AppendSeparator();
-
-    // Add the automatic rotation sub-menu
-    append_menu_item(&m_sla_object_menu, wxID_ANY, _L("Optimize orientation"), _L("Optimize the rotation of the object for better print results."),
-        [](wxCommandEvent&) { plater()->optimize_rotation(); });
 }
 
 void MenuFactory::create_part_menu()
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp
index c5060a88e..6b6905e4d 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.cpp
@@ -7,9 +7,10 @@
 
 #include "slic3r/GUI/GUI_App.hpp"
 #include "slic3r/GUI/GUI.hpp"
+#include "slic3r/GUI/Plater.hpp"
 #include "libslic3r/PresetBundle.hpp"
 
-#include "libslic3r/SLA/Rotfinder.hpp"
+#include "slic3r/GUI/Jobs/RotoptimizeJob.hpp"
 
 namespace Slic3r {
 namespace GUI {
@@ -204,6 +205,30 @@ void GLGizmoRotate3D::on_render_input_window(float x, float y, float bottom_limi
 {
     if (wxGetApp().preset_bundle->printers.get_edited_preset().printer_technology() != ptSLA)
         return;
+
+    RotoptimzeWindow popup{m_imgui, m_rotoptimizewin_state, {x, y, bottom_limit}};
+}
+
+void GLGizmoRotate3D::load_rotoptimize_state()
+{
+    std::string accuracy_str =
+        wxGetApp().app_config->get("sla_auto_rotate", "accuracy");
+
+    std::string method_str =
+        wxGetApp().app_config->get("sla_auto_rotate", "method_id");
+
+    if (!accuracy_str.empty()) {
+        float accuracy = std::stof(accuracy_str);
+        accuracy = std::max(0.f, std::min(accuracy, 1.f));
+
+        m_rotoptimizewin_state.accuracy = accuracy;
+    }
+
+    if (!method_str.empty()) {
+        int method_id = std::stoi(method_str);
+        if (method_id < int(RotoptimizeJob::get_methods_count()))
+            m_rotoptimizewin_state.method_id = method_id;
+    }
 }
 
 void GLGizmoRotate::render_circle() const
@@ -436,6 +461,8 @@ GLGizmoRotate3D::GLGizmoRotate3D(GLCanvas3D& parent, const std::string& icon_fil
     {
         m_gizmos[i].set_group_id(i);
     }
+
+    load_rotoptimize_state();
 }
 
 bool GLGizmoRotate3D::on_init()
@@ -492,5 +519,51 @@ void GLGizmoRotate3D::on_render() const
         m_gizmos[Z].render();
 }
 
+const char * GLGizmoRotate3D::RotoptimzeWindow::options[RotoptimizeJob::get_methods_count()];
+bool GLGizmoRotate3D::RotoptimzeWindow::options_valid = false;
+
+GLGizmoRotate3D::RotoptimzeWindow::RotoptimzeWindow(ImGuiWrapper *   imgui,
+                                                    State &          state,
+                                                    const Alignment &alignment)
+    : m_imgui{imgui}
+{
+    imgui->begin(_L("Optimize orientation"), 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);
+
+    ImGui::PushItemWidth(200.f);
+
+    size_t methods_cnt = RotoptimizeJob::get_methods_count();
+    if (!options_valid) {
+        for (size_t i = 0; i < methods_cnt; ++i)
+            options[i] = RotoptimizeJob::get_method_names()[i].c_str();
+
+        options_valid = true;
+    }
+
+    int citem = state.method_id;
+    if (ImGui::Combo(_L("Choose goal").c_str(), &citem, options, methods_cnt) ) {
+        state.method_id = citem;
+        wxGetApp().app_config->set("sla_auto_rotate", "method_id", std::to_string(state.method_id));
+    }
+
+    ImGui::Separator();
+
+    if ( imgui->button(_L("Optimize")) ) {
+        wxGetApp().plater()->optimize_rotation();
+    }
+}
+
+GLGizmoRotate3D::RotoptimzeWindow::~RotoptimzeWindow()
+{
+    m_imgui->end();
+}
+
 } // namespace GUI
 } // namespace Slic3r
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp
index 126c97b1d..c18d0eefd 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoRotate.hpp
@@ -2,6 +2,7 @@
 #define slic3r_GLGizmoRotate_hpp_
 
 #include "GLGizmoBase.hpp"
+#include "../Jobs/RotoptimizeJob.hpp"
 
 
 namespace Slic3r {
@@ -136,6 +137,39 @@ protected:
     }
 
     void on_render_input_window(float x, float y, float bottom_limit) override;
+
+private:
+
+    class RotoptimzeWindow {
+        ImGuiWrapper *m_imgui = nullptr;
+
+        static const char * options [];
+        static bool options_valid;
+
+    public:
+
+        struct State {
+            float  accuracy  = 1.f;
+            int    method_id = 0;
+        };
+
+        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 = {};
+
+    void load_rotoptimize_state();
 };
 
 } // namespace GUI
diff --git a/src/slic3r/GUI/Jobs/ArrangeJob.cpp b/src/slic3r/GUI/Jobs/ArrangeJob.cpp
index 3f1207b47..391c1fe28 100644
--- a/src/slic3r/GUI/Jobs/ArrangeJob.cpp
+++ b/src/slic3r/GUI/Jobs/ArrangeJob.cpp
@@ -78,7 +78,7 @@ void ArrangeJob::prepare_all() {
     for (ModelObject *obj: m_plater->model().objects)
         for (ModelInstance *mi : obj->instances) {
             ArrangePolygons & cont = mi->printable ? m_selected : m_unprintable;
-            cont.emplace_back(get_arrange_poly(PtrWrapper{mi}, m_plater));
+            cont.emplace_back(get_arrange_poly(mi, m_plater));
         }
 
     if (auto wti = get_wipe_tower_arrangepoly(*m_plater))
@@ -111,7 +111,7 @@ void ArrangeJob::prepare_selected() {
         
         for (size_t i = 0; i < inst_sel.size(); ++i) {
             ArrangePolygon &&ap =
-                get_arrange_poly(PtrWrapper{mo->instances[i]}, m_plater);
+                get_arrange_poly(mo->instances[i], m_plater);
 
             ArrangePolygons &cont = mo->instances[i]->printable ?
                         (inst_sel[i] ? m_selected :
@@ -161,12 +161,7 @@ void ArrangeJob::process()
 {
     static const auto arrangestr = _(L("Arranging"));
 
-    const GLCanvas3D::ArrangeSettings &settings =
-        static_cast<const GLCanvas3D*>(m_plater->canvas3D())->get_arrange_settings();
-
-    arrangement::ArrangeParams params;
-    params.allow_rotations  = settings.enable_rotation;
-    params.min_obj_distance = scaled(settings.distance);
+    arrangement::ArrangeParams params = get_arrange_params(m_plater);
 
     auto count = unsigned(m_selected.size() + m_unprintable.size());
     Points bedpts = get_bed_shape(*m_plater->config());
@@ -235,4 +230,23 @@ double bed_stride(const Plater *plater) {
     return scaled<double>((1. + LOGICAL_BED_GAP) * bedwidth);
 }
 
+template<>
+arrangement::ArrangePolygon get_arrange_poly(ModelInstance *inst,
+                                             const Plater * plater)
+{
+    return get_arrange_poly(PtrWrapper{inst}, plater);
+}
+
+arrangement::ArrangeParams get_arrange_params(Plater *p)
+{
+    const GLCanvas3D::ArrangeSettings &settings =
+        static_cast<const GLCanvas3D*>(p->canvas3D())->get_arrange_settings();
+
+    arrangement::ArrangeParams params;
+    params.allow_rotations  = settings.enable_rotation;
+    params.min_obj_distance = scaled(settings.distance);
+
+    return params;
+}
+
 }} // namespace Slic3r::GUI
diff --git a/src/slic3r/GUI/Jobs/ArrangeJob.hpp b/src/slic3r/GUI/Jobs/ArrangeJob.hpp
index 38aafb52c..7fa4b927e 100644
--- a/src/slic3r/GUI/Jobs/ArrangeJob.hpp
+++ b/src/slic3r/GUI/Jobs/ArrangeJob.hpp
@@ -4,7 +4,11 @@
 #include "PlaterJob.hpp"
 #include "libslic3r/Arrange.hpp"
 
-namespace Slic3r { namespace GUI {
+namespace Slic3r {
+
+class ModelInstance;
+
+namespace GUI {
 
 class ArrangeJob : public PlaterJob
 {
@@ -89,6 +93,11 @@ arrangement::ArrangePolygon get_arrange_poly(T obj, const Plater *plater)
     return ap;
 }
 
+template<>
+arrangement::ArrangePolygon get_arrange_poly(ModelInstance *inst,
+                                             const Plater * plater);
+
+arrangement::ArrangeParams get_arrange_params(Plater *p);
 
 }} // namespace Slic3r::GUI
 
diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp
index 978ccf8fc..a670affef 100644
--- a/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp
+++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.cpp
@@ -7,47 +7,84 @@
 #include "libslic3r/SLAPrint.hpp"
 
 #include "slic3r/GUI/Plater.hpp"
+#include "libslic3r/PresetBundle.hpp"
+
+#include "slic3r/GUI/GUI_App.hpp"
+#include "libslic3r/AppConfig.hpp"
 
 namespace Slic3r { namespace GUI {
 
+void RotoptimizeJob::prepare()
+{
+    std::string accuracy_str =
+        wxGetApp().app_config->get("sla_auto_rotate", "accuracy");
+
+    std::string method_str =
+        wxGetApp().app_config->get("sla_auto_rotate", "method_id");
+
+    if (!accuracy_str.empty())
+        m_accuracy = std::stof(accuracy_str);
+
+    if (!method_str.empty())
+        m_method_id = std::stoi(method_str);
+
+    m_accuracy = std::max(0.f, std::min(m_accuracy, 1.f));
+    m_method_id = std::max(size_t(0), std::min(get_methods_count() - 1, m_method_id));
+
+    m_default_print_cfg = wxGetApp().preset_bundle->full_config();
+
+    const auto &sel = m_plater->get_selection().get_content();
+
+    m_selected_object_ids.clear();
+    m_selected_object_ids.reserve(sel.size());
+    for (auto &[obj_idx, ignore] : sel)
+        m_selected_object_ids.emplace_back(obj_idx);
+}
+
 void RotoptimizeJob::process()
 {
-    int obj_idx = m_plater->get_selected_object_idx();
-    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)];
+    int prev_status = 0;
+    auto params =
+        sla::RotOptimizeParams{}
+            .accuracy(m_accuracy)
+            .print_config(&m_default_print_cfg)
+            .statucb([this, &prev_status](int s)
+        {
+            if (s > 0 && s < 100)
+                update_status(prev_status + s / m_selected_object_ids.size(),
+                              _(L("Searching for optimal orientation")));
 
-    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")));
-        },
-        [this] () { return was_canceled(); });
+            return !was_canceled();
+        });
 
 
-    double mindist = 6.0; // FIXME
+    for (ObjRot &objrot : m_selected_object_ids) {
+        ModelObject *o = m_plater->model().objects[size_t(objrot.idx)];
+        if (!o) continue;
+
+        if (Methods[m_method_id].findfn)
+            objrot.rot = Methods[m_method_id].findfn(*o, params);
+
+        prev_status += 100 / m_selected_object_ids.size();
+
+        if (was_canceled()) break;
+    }
+
+    update_status(100, was_canceled() ? _(L("Orientation search canceled.")) :
+                                        _(L("Orientation found.")));
+}
+
+void RotoptimizeJob::finalize()
+{
+    if (was_canceled()) return;
+
+    for (const ObjRot &objrot : m_selected_object_ids) {
+        ModelObject *o = m_plater->model().objects[size_t(objrot.idx)];
+        if (!o) continue;
 
-    if (!was_canceled()) {
         for(ModelInstance * oi : o->instances) {
-            oi->set_rotation({r[X], r[Y], 0.});
+            if (objrot.rot)
+                oi->set_rotation({objrot.rot->x(), objrot.rot->y(), 0.});
 
             auto    trmatrix = oi->get_transformation().get_matrix();
             Polygon trchull  = o->convex_hull_2d(trmatrix);
@@ -63,19 +100,13 @@ void RotoptimizeJob::process()
             oi->set_rotation(rt);
         }
 
-        m_plater->find_new_position(o->instances, scaled(mindist));
-
         // Correct the z offset of the object which was corrupted be
         // the rotation
         o->ensure_on_bed();
+
+//        m_plater->find_new_position(o->instances);
     }
 
-    update_status(100, was_canceled() ? _(L("Orientation search canceled.")) :
-                                        _(L("Orientation found.")));
-}
-
-void RotoptimizeJob::finalize()
-{
     if (!was_canceled())
         m_plater->update();
     
diff --git a/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp b/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp
index 06688b52d..3144f3c3e 100644
--- a/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp
+++ b/src/slic3r/GUI/Jobs/RotoptimizeJob.hpp
@@ -3,17 +3,70 @@
 
 #include "PlaterJob.hpp"
 
-namespace Slic3r { namespace GUI {
+#include "libslic3r/SLA/Rotfinder.hpp"
+#include "libslic3r/PrintConfig.hpp"
+
+namespace Slic3r {
+
+namespace GUI {
 
 class RotoptimizeJob : public PlaterJob
 {
+    using FindFn = std::function<Vec2d(const ModelObject &           mo,
+                                       const sla::RotOptimizeParams &params)>;
+
+    struct FindMethod { std::string name; FindFn findfn; };
+
+    static inline const FindMethod Methods[] = {
+        { L("Best surface quality"), sla::find_best_misalignment_rotation },
+        { L("Least supports"), sla::find_least_supports_rotation },
+        // Just a min area bounding box that is done for all methods anyway.
+        { L("Z axis only"), nullptr }
+    };
+
+    size_t m_method_id = 0;
+    float  m_accuracy  = 0.75;
+
+    DynamicPrintConfig m_default_print_cfg;
+
+    struct ObjRot
+    {
+        size_t               idx;
+        std::optional<Vec2d> rot;
+        ObjRot(size_t id): idx{id}, rot{} {}
+    };
+
+    std::vector<ObjRot> m_selected_object_ids;
+
+protected:
+
+    void prepare() override;
+
 public:
+
     RotoptimizeJob(std::shared_ptr<ProgressIndicator> pri, Plater *plater)
         : PlaterJob{std::move(pri), plater}
     {}
     
     void process() override;
     void finalize() override;
+
+    static constexpr size_t get_methods_count() { return std::size(Methods); }
+    static const auto & get_method_names()
+    {
+        static bool m_method_names_valid = false;
+        static std::array<std::string, std::size(Methods)> m_method_names;
+
+        if (!m_method_names_valid) {
+
+            for (size_t i = 0; i < std::size(Methods); ++i)
+                m_method_names[i] = _utf8(Methods[i].name);
+
+            m_method_names_valid = true;
+        }
+
+        return m_method_names;
+    }
 };
 
 }} // namespace Slic3r::GUI
diff --git a/src/slic3r/GUI/Plater.cpp b/src/slic3r/GUI/Plater.cpp
index a640d6e28..e0a5031c2 100644
--- a/src/slic3r/GUI/Plater.cpp
+++ b/src/slic3r/GUI/Plater.cpp
@@ -2669,32 +2669,36 @@ void Plater::priv::mirror(Axis axis)
     view3D->mirror_selection(axis);
 }
 
-void Plater::find_new_position(const ModelInstancePtrs &instances,
-                                     coord_t min_d)
+void Plater::find_new_position(const ModelInstancePtrs &instances)
 {
     arrangement::ArrangePolygons movable, fixed;
+    arrangement::ArrangeParams arr_params = get_arrange_params(this);
     
     for (const ModelObject *mo : p->model.objects)
-        for (const ModelInstance *inst : mo->instances) {
+        for (ModelInstance *inst : mo->instances) {
             auto it = std::find(instances.begin(), instances.end(), inst);
-            auto arrpoly = inst->get_arrange_polygon();
+            auto arrpoly = get_arrange_poly(inst, this);
 
             if (it == instances.end())
                 fixed.emplace_back(std::move(arrpoly));
-            else
+            else {
+                arrpoly.setter = [it](const arrangement::ArrangePolygon &p) {
+                    if (p.is_arranged() && p.bed_idx == 0) {
+                        Vec2d t = p.translation.cast<double>();
+                        (*it)->apply_arrange_result(t, p.rotation);
+                    }
+                };
                 movable.emplace_back(std::move(arrpoly));
+            }
         }
     
     if (auto wt = get_wipe_tower_arrangepoly(*this))
         fixed.emplace_back(*wt);
     
-    arrangement::arrange(movable, fixed, get_bed_shape(*config()),
-                         arrangement::ArrangeParams{min_d});
+    arrangement::arrange(movable, fixed, get_bed_shape(*config()), arr_params);
 
-    for (size_t i = 0; i < instances.size(); ++i)
-        if (movable[i].bed_idx == 0)
-            instances[i]->apply_arrange_result(movable[i].translation.cast<double>(),
-                                               movable[i].rotation);
+    for (auto & m : movable)
+        m.apply();
 }
 
 void Plater::priv::split_object()
diff --git a/src/slic3r/GUI/Plater.hpp b/src/slic3r/GUI/Plater.hpp
index ff81dad26..16590ed9b 100644
--- a/src/slic3r/GUI/Plater.hpp
+++ b/src/slic3r/GUI/Plater.hpp
@@ -258,7 +258,7 @@ public:
     BoundingBoxf bed_shape_bb() const;
     
     void arrange();
-    void find_new_position(const ModelInstancePtrs  &instances, coord_t min_d);
+    void find_new_position(const ModelInstancePtrs  &instances);
 
     void set_current_canvas_as_dirty();
     void unbind_canvas_event_handlers();
diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp
index bdd5731dc..1f98463cc 100644
--- a/tests/sla_print/sla_print_tests.cpp
+++ b/tests/sla_print/sla_print_tests.cpp
@@ -1,6 +1,7 @@
 #include <unordered_set>
 #include <unordered_map>
 #include <random>
+#include <numeric>
 #include <cstdint>
 
 #include "sla_test_utils.hpp"
@@ -248,7 +249,7 @@ TEST_CASE("Test concurrency")
 
     double ref = std::accumulate(vals.begin(), vals.end(), 0.);
 
-    double s = sla::ccr_par::reduce(vals.begin(), vals.end(), 0., std::plus<double>{});
+    double s = execution::accumulate(ex_tbb, vals.begin(), vals.end(), 0.);
 
     REQUIRE(s == Approx(ref));
 }