diff --git a/src/admesh/stl.h b/src/admesh/stl.h
index 9224b0459..e0f2865f0 100644
--- a/src/admesh/stl.h
+++ b/src/admesh/stl.h
@@ -255,18 +255,24 @@ extern void its_transform(indexed_triangle_set &its, T *trafo3x4)
 }
 
 template<typename T>
-inline void its_transform(indexed_triangle_set &its, const Eigen::Transform<T, 3, Eigen::Affine, Eigen::DontAlign>& t)
+inline void its_transform(indexed_triangle_set &its, const Eigen::Transform<T, 3, Eigen::Affine, Eigen::DontAlign>& t, bool fix_left_handed = false)
 {
 	//const Eigen::Matrix<double, 3, 3, Eigen::DontAlign> r = t.matrix().template block<3, 3>(0, 0);
 	for (stl_vertex &v : its.vertices)
 		v = (t * v.template cast<T>()).template cast<float>().eval();
+  if (fix_left_handed && t.matrix().block(0, 0, 3, 3).determinant() < 0.)
+    for (stl_triangle_vertex_indices &i : its.indices)
+      std::swap(i[0], i[1]);
 }
 
 template<typename T>
-inline void its_transform(indexed_triangle_set &its, const Eigen::Matrix<T, 3, 3, Eigen::DontAlign>& m)
+inline void its_transform(indexed_triangle_set &its, const Eigen::Matrix<T, 3, 3, Eigen::DontAlign>& m, bool fix_left_handed = false)
 {
-	for (stl_vertex &v : its.vertices)
+  for (stl_vertex &v : its.vertices)
 		v = (m * v.template cast<T>()).template cast<float>().eval();
+  if (fix_left_handed && m.determinant() < 0.)
+    for (stl_triangle_vertex_indices &i : its.indices)
+      std::swap(i[0], i[1]);
 }
 
 extern void its_rotate_x(indexed_triangle_set &its, float angle);
diff --git a/src/libslic3r/AABBTreeIndirect.hpp b/src/libslic3r/AABBTreeIndirect.hpp
index 17d918aeb..964133faa 100644
--- a/src/libslic3r/AABBTreeIndirect.hpp
+++ b/src/libslic3r/AABBTreeIndirect.hpp
@@ -283,7 +283,7 @@ namespace detail {
 
 	template<typename V, typename W>
     std::enable_if_t<std::is_same<typename V::Scalar, double>::value && std::is_same<typename W::Scalar, double>::value, bool>
-	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W v2, double &t, double &u, double &v) {
+	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) {
         return intersect_triangle1(const_cast<double*>(origin.data()), const_cast<double*>(dir.data()),
                                    const_cast<double*>(v0.data()), const_cast<double*>(v1.data()), const_cast<double*>(v2.data()),
                                    &t, &u, &v);
@@ -291,7 +291,7 @@ namespace detail {
 
 	template<typename V, typename W>
     std::enable_if_t<std::is_same<typename V::Scalar, double>::value && !std::is_same<typename W::Scalar, double>::value, bool>
-	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W v2, double &t, double &u, double &v) {
+	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) {
         using Vector = Eigen::Matrix<double, 3, 1>;
         Vector w0 = v0.template cast<double>();
         Vector w1 = v1.template cast<double>();
@@ -302,7 +302,7 @@ namespace detail {
 
 	template<typename V, typename W>
     std::enable_if_t<! std::is_same<typename V::Scalar, double>::value && std::is_same<typename W::Scalar, double>::value, bool>
-	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W v2, double &t, double &u, double &v) {
+	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) {
         using Vector = Eigen::Matrix<double, 3, 1>;
         Vector o  = origin.template cast<double>();
         Vector d  = dir.template cast<double>();
@@ -311,7 +311,7 @@ namespace detail {
 
 	template<typename V, typename W>
     std::enable_if_t<! std::is_same<typename V::Scalar, double>::value && ! std::is_same<typename W::Scalar, double>::value, bool>
-	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W v2, double &t, double &u, double &v) {
+	intersect_triangle(const V &origin, const V &dir, const W &v0, const W &v1, const W &v2, double &t, double &u, double &v) {
         using Vector = Eigen::Matrix<double, 3, 1>;
         Vector o  = origin.template cast<double>();
         Vector d  = dir.template cast<double>();
diff --git a/src/libslic3r/BoundingBox.cpp b/src/libslic3r/BoundingBox.cpp
index e3f939509..eb4e042a0 100644
--- a/src/libslic3r/BoundingBox.cpp
+++ b/src/libslic3r/BoundingBox.cpp
@@ -75,6 +75,7 @@ BoundingBoxBase<PointClass>::merge(const PointClass &point)
     }
 }
 template void BoundingBoxBase<Point>::merge(const Point &point);
+template void BoundingBoxBase<Vec2f>::merge(const Vec2f &point);
 template void BoundingBoxBase<Vec2d>::merge(const Vec2d &point);
 
 template <class PointClass> void
@@ -101,6 +102,7 @@ BoundingBoxBase<PointClass>::merge(const BoundingBoxBase<PointClass> &bb)
     }
 }
 template void BoundingBoxBase<Point>::merge(const BoundingBoxBase<Point> &bb);
+template void BoundingBoxBase<Vec2f>::merge(const BoundingBoxBase<Vec2f> &bb);
 template void BoundingBoxBase<Vec2d>::merge(const BoundingBoxBase<Vec2d> &bb);
 
 template <class PointClass> void
@@ -115,6 +117,7 @@ BoundingBox3Base<PointClass>::merge(const PointClass &point)
         this->defined = true;
     }
 }
+template void BoundingBox3Base<Vec3f>::merge(const Vec3f &point);
 template void BoundingBox3Base<Vec3d>::merge(const Vec3d &point);
 
 template <class PointClass> void
@@ -147,6 +150,7 @@ BoundingBoxBase<PointClass>::size() const
     return PointClass(this->max(0) - this->min(0), this->max(1) - this->min(1));
 }
 template Point BoundingBoxBase<Point>::size() const;
+template Vec2f BoundingBoxBase<Vec2f>::size() const;
 template Vec2d BoundingBoxBase<Vec2d>::size() const;
 
 template <class PointClass> PointClass
@@ -154,6 +158,7 @@ BoundingBox3Base<PointClass>::size() const
 {
     return PointClass(this->max(0) - this->min(0), this->max(1) - this->min(1), this->max(2) - this->min(2));
 }
+template Vec3f BoundingBox3Base<Vec3f>::size() const;
 template Vec3d BoundingBox3Base<Vec3d>::size() const;
 
 template <class PointClass> double BoundingBoxBase<PointClass>::radius() const
@@ -200,6 +205,7 @@ BoundingBoxBase<PointClass>::center() const
     return (this->min + this->max) / 2;
 }
 template Point BoundingBoxBase<Point>::center() const;
+template Vec2f BoundingBoxBase<Vec2f>::center() const;
 template Vec2d BoundingBoxBase<Vec2d>::center() const;
 
 template <class PointClass> PointClass
@@ -207,6 +213,7 @@ BoundingBox3Base<PointClass>::center() const
 {
     return (this->min + this->max) / 2;
 }
+template Vec3f BoundingBox3Base<Vec3f>::center() const;
 template Vec3d BoundingBox3Base<Vec3d>::center() const;
 
 template <class PointClass> coordf_t
@@ -215,6 +222,7 @@ BoundingBox3Base<PointClass>::max_size() const
     PointClass s = size();
     return std::max(s(0), std::max(s(1), s(2)));
 }
+template coordf_t BoundingBox3Base<Vec3f>::max_size() const;
 template coordf_t BoundingBox3Base<Vec3d>::max_size() const;
 
 // Align a coordinate to a grid. The coordinate may be negative,
diff --git a/src/libslic3r/BoundingBox.hpp b/src/libslic3r/BoundingBox.hpp
index 71746f32e..065476cb2 100644
--- a/src/libslic3r/BoundingBox.hpp
+++ b/src/libslic3r/BoundingBox.hpp
@@ -19,6 +19,8 @@ public:
     BoundingBoxBase() : min(PointClass::Zero()), max(PointClass::Zero()), defined(false) {}
     BoundingBoxBase(const PointClass &pmin, const PointClass &pmax) : 
         min(pmin), max(pmax), defined(pmin(0) < pmax(0) && pmin(1) < pmax(1)) {}
+    BoundingBoxBase(const PointClass &p1, const PointClass &p2, const PointClass &p3) :
+        min(p1), max(p1), defined(false) { merge(p2); merge(p3); }
     BoundingBoxBase(const std::vector<PointClass>& points) : min(PointClass::Zero()), max(PointClass::Zero())
     {
         if (points.empty()) {
@@ -66,6 +68,8 @@ public:
     BoundingBox3Base(const PointClass &pmin, const PointClass &pmax) : 
         BoundingBoxBase<PointClass>(pmin, pmax) 
         { if (pmin(2) >= pmax(2)) BoundingBoxBase<PointClass>::defined = false; }
+    BoundingBox3Base(const PointClass &p1, const PointClass &p2, const PointClass &p3) :
+        BoundingBoxBase<PointClass>(p1, p1) { merge(p2); merge(p3); }
     BoundingBox3Base(const std::vector<PointClass>& points)
     {
         if (points.empty())
@@ -110,24 +114,32 @@ extern template void     BoundingBoxBase<Vec3d>::scale(double factor);
 extern template void     BoundingBoxBase<Point>::offset(coordf_t delta);
 extern template void     BoundingBoxBase<Vec2d>::offset(coordf_t delta);
 extern template void     BoundingBoxBase<Point>::merge(const Point &point);
+extern template void     BoundingBoxBase<Vec2f>::merge(const Vec2f &point);
 extern template void     BoundingBoxBase<Vec2d>::merge(const Vec2d &point);
 extern template void     BoundingBoxBase<Point>::merge(const Points &points);
 extern template void     BoundingBoxBase<Vec2d>::merge(const Pointfs &points);
 extern template void     BoundingBoxBase<Point>::merge(const BoundingBoxBase<Point> &bb);
+extern template void     BoundingBoxBase<Vec2f>::merge(const BoundingBoxBase<Vec2f> &bb);
 extern template void     BoundingBoxBase<Vec2d>::merge(const BoundingBoxBase<Vec2d> &bb);
 extern template Point    BoundingBoxBase<Point>::size() const;
+extern template Vec2f    BoundingBoxBase<Vec2f>::size() const;
 extern template Vec2d    BoundingBoxBase<Vec2d>::size() const;
 extern template double   BoundingBoxBase<Point>::radius() const;
 extern template double   BoundingBoxBase<Vec2d>::radius() const;
 extern template Point    BoundingBoxBase<Point>::center() const;
+extern template Vec2f    BoundingBoxBase<Vec2f>::center() const;
 extern template Vec2d    BoundingBoxBase<Vec2d>::center() const;
+extern template void     BoundingBox3Base<Vec3f>::merge(const Vec3f &point);
 extern template void     BoundingBox3Base<Vec3d>::merge(const Vec3d &point);
 extern template void     BoundingBox3Base<Vec3d>::merge(const Pointf3s &points);
 extern template void     BoundingBox3Base<Vec3d>::merge(const BoundingBox3Base<Vec3d> &bb);
+extern template Vec3f    BoundingBox3Base<Vec3f>::size() const;
 extern template Vec3d    BoundingBox3Base<Vec3d>::size() const;
 extern template double   BoundingBox3Base<Vec3d>::radius() const;
 extern template void     BoundingBox3Base<Vec3d>::offset(coordf_t delta);
+extern template Vec3f    BoundingBox3Base<Vec3f>::center() const;
 extern template Vec3d    BoundingBox3Base<Vec3d>::center() const;
+extern template coordf_t BoundingBox3Base<Vec3f>::max_size() const;
 extern template coordf_t BoundingBox3Base<Vec3d>::max_size() const;
 
 class BoundingBox : public BoundingBoxBase<Point>
diff --git a/src/libslic3r/Fill/Fill.cpp b/src/libslic3r/Fill/Fill.cpp
index d68bc7afb..70792b823 100644
--- a/src/libslic3r/Fill/Fill.cpp
+++ b/src/libslic3r/Fill/Fill.cpp
@@ -345,8 +345,7 @@ void Layer::make_fills(FillAdaptive_Internal::Octree* adaptive_fill_octree, Fill
         f->layer_id = this->id();
         f->z 		= this->print_z;
         f->angle 	= surface_fill.params.angle;
-        f->adapt_fill_octree = adaptive_fill_octree;
-        f->support_fill_octree = support_fill_octree;
+        f->adapt_fill_octree = (surface_fill.params.pattern == ipSupportCubic) ? support_fill_octree : adaptive_fill_octree;
 
         // calculate flow spacing for infill pattern generation
         bool using_internal_flow = ! surface_fill.surface.is_solid() && ! surface_fill.params.flow.bridge;
diff --git a/src/libslic3r/Fill/FillAdaptive.cpp b/src/libslic3r/Fill/FillAdaptive.cpp
index 3b9212230..6c5d7c0c4 100644
--- a/src/libslic3r/Fill/FillAdaptive.cpp
+++ b/src/libslic3r/Fill/FillAdaptive.cpp
@@ -2,16 +2,211 @@
 #include "../ExPolygon.hpp"
 #include "../Surface.hpp"
 #include "../Geometry.hpp"
-#include "../AABBTreeIndirect.hpp"
 #include "../Layer.hpp"
 #include "../Print.hpp"
 #include "../ShortestPath.hpp"
 
 #include "FillAdaptive.hpp"
 
+// for indexed_triangle_set
+#include <admesh/stl.h>
+
+#include <cstdlib>
+#include <cmath>
+
+// Boost pool: Don't use mutexes to synchronize memory allocation.
+#define BOOST_POOL_NO_MT
+#include <boost/pool/object_pool.hpp>
+
 namespace Slic3r {
 
-std::pair<double, double> adaptive_fill_line_spacing(const PrintObject &print_object)
+// Derived from https://github.com/juj/MathGeoLib/blob/master/src/Geometry/Triangle.cpp
+// The AABB-Triangle test implementation is based on the pseudo-code in
+// Christer Ericson's Real-Time Collision Detection, pp. 169-172. It is
+// practically a standard SAT test.
+//
+// Original MathGeoLib benchmark:
+//    Best: 17.282 nsecs / 46.496 ticks, Avg: 17.804 nsecs, Worst: 18.434 nsecs
+//
+//FIXME Vojtech: The MathGeoLib contains a vectorized implementation.
+template<typename Vector> 
+bool triangle_AABB_intersects(const Vector &a, const Vector &b, const Vector &c, const BoundingBoxBase<Vector> &aabb)
+{
+    using Scalar = typename Vector::Scalar;
+
+    Vector tMin = a.cwiseMin(b.cwiseMin(c));
+    Vector tMax = a.cwiseMax(b.cwiseMax(c));
+
+    if (tMin.x() >= aabb.max.x() || tMax.x() <= aabb.min.x()
+        || tMin.y() >= aabb.max.y() || tMax.y() <= aabb.min.y()
+        || tMin.z() >= aabb.max.z() || tMax.z() <= aabb.min.z())
+        return false;
+
+    Vector center = (aabb.min + aabb.max) * 0.5f;
+    Vector h = aabb.max - center;
+
+    const Vector t[3] { b-a, c-a, c-b };
+
+    Vector ac = a - center;
+
+    Vector n = t[0].cross(t[1]);
+    Scalar s = n.dot(ac);
+    Scalar r = std::abs(h.dot(n.cwiseAbs()));
+    if (abs(s) >= r)
+        return false;
+
+    const Vector at[3] = { t[0].cwiseAbs(), t[1].cwiseAbs(), t[2].cwiseAbs() };
+
+    Vector bc = b - center;
+    Vector cc = c - center;
+
+    // SAT test all cross-axes.
+    // The following is a fully unrolled loop of this code, stored here for reference:
+    /*
+    Scalar d1, d2, a1, a2;
+    const Vector e[3] = { DIR_VEC(1, 0, 0), DIR_VEC(0, 1, 0), DIR_VEC(0, 0, 1) };
+    for(int i = 0; i < 3; ++i)
+        for(int j = 0; j < 3; ++j)
+        {
+            Vector axis = Cross(e[i], t[j]);
+            ProjectToAxis(axis, d1, d2);
+            aabb.ProjectToAxis(axis, a1, a2);
+            if (d2 <= a1 || d1 >= a2) return false;
+        }
+    */
+
+    // eX <cross> t[0]
+    Scalar d1 = t[0].y() * ac.z() - t[0].z() * ac.y();
+    Scalar d2 = t[0].y() * cc.z() - t[0].z() * cc.y();
+    Scalar tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[0].z() + h.z() * at[0].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eX <cross> t[1]
+    d1 = t[1].y() * ac.z() - t[1].z() * ac.y();
+    d2 = t[1].y() * bc.z() - t[1].z() * bc.y();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[1].z() + h.z() * at[1].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eX <cross> t[2]
+    d1 = t[2].y() * ac.z() - t[2].z() * ac.y();
+    d2 = t[2].y() * bc.z() - t[2].z() * bc.y();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[2].z() + h.z() * at[2].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eY <cross> t[0]
+    d1 = t[0].z() * ac.x() - t[0].x() * ac.z();
+    d2 = t[0].z() * cc.x() - t[0].x() * cc.z();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.x() * at[0].z() + h.z() * at[0].x());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eY <cross> t[1]
+    d1 = t[1].z() * ac.x() - t[1].x() * ac.z();
+    d2 = t[1].z() * bc.x() - t[1].x() * bc.z();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.x() * at[1].z() + h.z() * at[1].x());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eY <cross> t[2]
+    d1 = t[2].z() * ac.x() - t[2].x() * ac.z();
+    d2 = t[2].z() * bc.x() - t[2].x() * bc.z();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.x() * at[2].z() + h.z() * at[2].x());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eZ <cross> t[0]
+    d1 = t[0].x() * ac.y() - t[0].y() * ac.x();
+    d2 = t[0].x() * cc.y() - t[0].y() * cc.x();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[0].x() + h.x() * at[0].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eZ <cross> t[1]
+    d1 = t[1].x() * ac.y() - t[1].y() * ac.x();
+    d2 = t[1].x() * bc.y() - t[1].y() * bc.x();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[1].x() + h.x() * at[1].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // eZ <cross> t[2]
+    d1 = t[2].x() * ac.y() - t[2].y() * ac.x();
+    d2 = t[2].x() * bc.y() - t[2].y() * bc.x();
+    tc = (d1 + d2) * 0.5f;
+    r = std::abs(h.y() * at[2].x() + h.x() * at[2].y());
+    if (r + std::abs(tc - d1) < std::abs(tc))
+        return false;
+
+    // No separating axis exists, the AABB and triangle intersect.
+    return true;
+}
+
+// Ordering of children cubes.
+static const std::array<Vec3d, 8> child_centers {
+    Vec3d(-1, -1, -1), Vec3d( 1, -1, -1), Vec3d(-1,  1, -1), Vec3d( 1,  1, -1),
+    Vec3d(-1, -1,  1), Vec3d( 1, -1,  1), Vec3d(-1,  1,  1), Vec3d( 1,  1,  1)
+};
+
+// Traversal order of octree children cells for three infill directions,
+// so that a single line will be discretized in a strictly monotonous order.
+static constexpr std::array<std::array<int, 8>, 3> child_traversal_order {
+    std::array<int, 8>{ 2, 3, 0, 1, 6, 7, 4, 5 },
+    std::array<int, 8>{ 4, 0, 6, 2, 5, 1, 7, 3 },
+    std::array<int, 8>{ 1, 5, 0, 4, 3, 7, 2, 6 },
+};
+
+namespace FillAdaptive_Internal
+{
+    struct Cube
+    {
+        Vec3d center;
+#ifndef NDEBUG
+        Vec3d center_octree;
+#endif // NDEBUG
+        std::array<Cube*, 8> children {}; // initialized to nullptrs
+        Cube(const Vec3d &center) : center(center) {}
+    };
+
+    struct CubeProperties
+    {
+        double edge_length;     // Lenght of edge of a cube
+        double height;          // Height of rotated cube (standing on the corner)
+        double diagonal_length; // Length of diagonal of a cube a face
+        double line_z_distance; // Defines maximal distance from a center of a cube on Z axis on which lines will be created
+        double line_xy_distance;// Defines maximal distance from a center of a cube on X and Y axis on which lines will be created
+    };
+
+    struct Octree
+    {
+        // Octree will allocate its Cubes from the pool. The pool only supports deletion of the complete pool,
+        // perfect for building up our octree.
+        boost::object_pool<Cube>    pool;
+        Cube*                       root_cube { nullptr };
+        Vec3d                       origin;
+        std::vector<CubeProperties> cubes_properties;
+
+        Octree(const Vec3d &origin, const std::vector<CubeProperties> &cubes_properties)
+            : root_cube(pool.construct(origin)), origin(origin), cubes_properties(cubes_properties) {}
+
+        void insert_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, Cube *current_cube, const BoundingBoxf3 &current_bbox, int depth);
+    };
+
+    void OctreeDeleter::operator()(Octree *p) {
+        delete p;
+    }
+}; // namespace FillAdaptive_Internal
+
+std::pair<double, double> FillAdaptive_Internal::adaptive_fill_line_spacing(const PrintObject &print_object)
 {
     // Output, spacing for icAdaptiveCubic and icSupportCubic
     double  adaptive_line_spacing = 0.;
@@ -90,431 +285,392 @@ std::pair<double, double> adaptive_fill_line_spacing(const PrintObject &print_ob
     return std::make_pair(adaptive_line_spacing, support_line_spacing);
 }
 
-void FillAdaptive::_fill_surface_single(const FillParams &             params,
-                                        unsigned int                   thickness_layers,
-                                        const std::pair<float, Point> &direction,
-                                        ExPolygon &                    expolygon,
-                                        Polylines &                    polylines_out)
+// Context used by generate_infill_lines() when recursively traversing an octree in a DDA fashion
+// (Digital Differential Analyzer).
+struct FillContext
 {
-    if(this->adapt_fill_octree != nullptr)
-        this->generate_infill(params, thickness_layers, direction, expolygon, polylines_out, this->adapt_fill_octree);
+    // The angles have to agree with child_traversal_order.
+    static constexpr double direction_angles[3] {
+        0.,
+        (2.0 * M_PI) / 3.0,
+        -(2.0 * M_PI) / 3.0
+    };
+
+    FillContext(const FillAdaptive_Internal::Octree &octree, double z_position, int direction_idx) :
+        origin_world(octree.origin),
+        cubes_properties(octree.cubes_properties),
+        z_position(z_position),
+        traversal_order(child_traversal_order[direction_idx]),
+        cos_a(cos(direction_angles[direction_idx])),
+        sin_a(sin(direction_angles[direction_idx]))
+    {
+        static constexpr auto unused = std::numeric_limits<coord_t>::max();
+        temp_lines.assign((1 << octree.cubes_properties.size()) - 1, Line(Point(unused, unused), Point(unused, unused)));
+    }
+
+    // Rotate the point, uses the same convention as Point::rotate().
+    Vec2d rotate(const Vec2d& v) { return Vec2d(this->cos_a * v.x() - this->sin_a * v.y(), this->sin_a * v.x() + this->cos_a * v.y()); }
+
+    // Center of the root cube in the Octree coordinate system.
+    const Vec3d                                                 origin_world;
+    const std::vector<FillAdaptive_Internal::CubeProperties>   &cubes_properties;
+    // Top of the current layer.
+    const double                                                z_position;
+    // Order of traversal for this line direction.
+    const std::array<int, 8>                                    traversal_order;
+    // Rotation of the generated line for this line direction.
+    const double                                                cos_a;
+    const double                                                sin_a;
+
+    // Linearized tree spanning a single Octree wall, used to connect lines spanning
+    // neighboring Octree cells. Unused lines have the Line::a::x set to infinity.
+    std::vector<Line>                                           temp_lines;
+    // Final output
+    std::vector<Line>                                           output_lines;
+};
+
+static constexpr double octree_rot[3] = { 5.0 * M_PI / 4.0, Geometry::deg2rad(215.264), M_PI / 6.0 };
+
+Eigen::Quaterniond FillAdaptive_Internal::adaptive_fill_octree_transform_to_world()
+{
+    return Eigen::AngleAxisd(octree_rot[2], Vec3d::UnitZ()) * Eigen::AngleAxisd(octree_rot[1], Vec3d::UnitY()) * Eigen::AngleAxisd(octree_rot[0], Vec3d::UnitX());
 }
 
-void FillAdaptive::generate_infill(const FillParams &             params,
-                                        unsigned int                   thickness_layers,
-                                        const std::pair<float, Point> &direction,
-                                        ExPolygon &                    expolygon,
-                                        Polylines &                    polylines_out,
-                                        FillAdaptive_Internal::Octree *octree)
+Eigen::Quaterniond FillAdaptive_Internal::adaptive_fill_octree_transform_to_octree()
 {
-    Vec3d rotation = Vec3d((5.0 * M_PI) / 4.0, Geometry::deg2rad(215.264), M_PI / 6.0);
-    Transform3d rotation_matrix = Geometry::assemble_transform(Vec3d::Zero(), rotation, Vec3d::Ones(), Vec3d::Ones());
+    return Eigen::AngleAxisd(- octree_rot[0], Vec3d::UnitX()) * Eigen::AngleAxisd(- octree_rot[1], Vec3d::UnitY()) * Eigen::AngleAxisd(- octree_rot[2], Vec3d::UnitZ());
+}
 
-    // Store grouped lines by its direction (multiple of 120°)
-    std::vector<Lines> infill_lines_dir(3);
-    this->generate_infill_lines(octree->root_cube.get(),
-                                this->z, octree->origin, rotation_matrix,
-                                infill_lines_dir, octree->cubes_properties,
-                                int(octree->cubes_properties.size()) - 1);
+#ifndef NDEBUG
+// Verify that the traversal order of the octree children matches the line direction,
+// therefore the infill line may get extended with O(1) time & space complexity.
+static bool verify_traversal_order(
+    FillContext                         &context,
+    const FillAdaptive_Internal::Cube   *cube,
+    int                                  depth,
+    const Vec2d                         &line_from,
+    const Vec2d                         &line_to)
+{
+    std::array<Vec3d, 8> c;
+    Eigen::Quaterniond to_world = FillAdaptive_Internal::adaptive_fill_octree_transform_to_world();
+    for (int i = 0; i < 8; ++i) {
+        int j = context.traversal_order[i];
+        Vec3d cntr = to_world * (cube->center_octree + (child_centers[j] * (context.cubes_properties[depth].edge_length / 4.)));
+        assert(!cube->children[j] || cube->children[j]->center.isApprox(cntr));
+        c[i] = cntr;
+    }
+    std::array<Vec3d, 10> dirs = {
+        c[1] - c[0], c[2] - c[0], c[3] - c[1], c[3] - c[2], c[3] - c[0],
+        c[5] - c[4], c[6] - c[4], c[7] - c[5], c[7] - c[6], c[7] - c[4]
+    };
+    assert(std::abs(dirs[4].z()) < 0.001);
+    assert(std::abs(dirs[9].z()) < 0.001);
+    assert(dirs[0].isApprox(dirs[3]));
+    assert(dirs[1].isApprox(dirs[2]));
+    assert(dirs[5].isApprox(dirs[8]));
+    assert(dirs[6].isApprox(dirs[7]));
+    Vec3d line_dir = Vec3d(line_to.x() - line_from.x(), line_to.y() - line_from.y(), 0.).normalized();
+    for (auto& dir : dirs) {
+        double d = dir.normalized().dot(line_dir);
+        assert(d > 0.7);
+    }
+    return true;
+}
+#endif // NDEBUG
+
+static void generate_infill_lines_recursive(
+    FillContext                         &context,
+    const FillAdaptive_Internal::Cube   *cube,
+    // Address of this wall in the octree,  used to address context.temp_lines.
+    int                                  address,
+    int                                  depth)
+{
+    assert(cube != nullptr);
+
+    const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties = context.cubes_properties;
+    const double z_diff     = context.z_position - cube->center.z();
+    const double z_diff_abs = std::abs(z_diff);
+
+    if (z_diff_abs > cubes_properties[depth].height / 2.)
+        return;
+
+    if (z_diff_abs < cubes_properties[depth].line_z_distance) {
+        // Discretize a single wall splitting the cube into two.
+        const double zdist = cubes_properties[depth].line_z_distance;
+        Vec2d from(
+            0.5 * cubes_properties[depth].diagonal_length * (zdist - z_diff_abs) / zdist,
+            cubes_properties[depth].line_xy_distance - (zdist + z_diff) / sqrt(2.));
+        Vec2d to(-from.x(), from.y());
+        from = context.rotate(from);
+        to   = context.rotate(to);
+        // Relative to cube center
+        Vec2d offset(cube->center.x() - context.origin_world.x(), cube->center.y() - context.origin_world.y());
+        from += offset;
+        to   += offset;
+        // Verify that the traversal order of the octree children matches the line direction,
+        // therefore the infill line may get extended with O(1) time & space complexity.
+        assert(verify_traversal_order(context, cube, depth, from, to));
+        // Either extend an existing line or start a new one.
+        Line &last_line = context.temp_lines[address];
+        Line  new_line(Point::new_scale(from), Point::new_scale(to));
+        if (last_line.a.x() == std::numeric_limits<coord_t>::max()) {
+            last_line.a = new_line.a;
+        } else if ((new_line.a - last_line.b).cwiseAbs().maxCoeff() > 300) { // SCALED_EPSILON is 100 and it is not enough) {
+            context.output_lines.emplace_back(last_line);
+            last_line.a = new_line.a;
+        }
+        last_line.b = new_line.b;
+    }
+
+    // left child index
+    address = address * 2 + 1;
+    -- depth;
+    size_t i = 0;
+    for (const int child_idx : context.traversal_order) {
+        const FillAdaptive_Internal::Cube *child = cube->children[child_idx];
+        if (child != nullptr)
+            generate_infill_lines_recursive(context, child, address, depth);
+        if (++ i == 4)
+            // right child index
+            ++ address;
+    }
+}
+
+#ifndef NDEBUG
+//    #define ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
+#endif
+
+#ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
+static void export_infill_lines_to_svg(const ExPolygon &expoly, const Polylines &polylines, const std::string &path)
+{
+    BoundingBox bbox = get_extents(expoly);
+    bbox.offset(scale_(3.));
+
+    ::Slic3r::SVG svg(path, bbox);
+    svg.draw(expoly);
+    svg.draw_outline(expoly, "green");
+    svg.draw(polylines, "red");
+    static constexpr double trim_length = scale_(0.4);
+    for (Polyline polyline : polylines) {
+        Vec2d a = polyline.points.front().cast<double>();
+        Vec2d d = polyline.points.back().cast<double>();
+        if (polyline.size() == 2) {
+            Vec2d v = d - a;
+            double l = v.norm();
+            if (l > 2. * trim_length) {
+                a += v * trim_length / l;
+                d -= v * trim_length / l;
+                polyline.points.front() = a.cast<coord_t>();
+                polyline.points.back() = d.cast<coord_t>();
+            } else
+                polyline.points.clear();
+        } else if (polyline.size() > 2) {
+            Vec2d b = polyline.points[1].cast<double>();
+            Vec2d c = polyline.points[polyline.points.size() - 2].cast<double>();
+            Vec2d v = b - a;
+            double l = v.norm();
+            if (l > trim_length) {
+                a += v * trim_length / l;
+                polyline.points.front() = a.cast<coord_t>();
+            } else
+                polyline.points.erase(polyline.points.begin());
+            v = d - c;
+            l = v.norm();
+            if (l > trim_length)
+                polyline.points.back() = (d - v * trim_length / l).cast<coord_t>();
+            else
+                polyline.points.pop_back();
+        }
+        svg.draw(polyline, "black");
+    }
+}
+#endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
+
+void FillAdaptive::_fill_surface_single(
+    const FillParams &             params,
+    unsigned int                   thickness_layers,
+    const std::pair<float, Point> &direction,
+    ExPolygon                     &expolygon,
+    Polylines                     &polylines_out)
+{
+    assert (this->adapt_fill_octree);
 
     Polylines all_polylines;
-    all_polylines.reserve(infill_lines_dir[0].size() * 3);
-    for (Lines &infill_lines : infill_lines_dir)
     {
-        for (const Line &line : infill_lines)
-        {
-            all_polylines.emplace_back(line.a, line.b);
+        // 3 contexts for three directions of infill lines
+        std::array<FillContext, 3> contexts { 
+            FillContext { *adapt_fill_octree, this->z, 0 },
+            FillContext { *adapt_fill_octree, this->z, 1 },
+            FillContext { *adapt_fill_octree, this->z, 2 }
+        };
+        // Generate the infill lines along the octree cells, merge touching lines of the same direction.
+        size_t num_lines = 0;
+        for (auto &context : contexts) {
+            generate_infill_lines_recursive(context, adapt_fill_octree->root_cube, 0, int(adapt_fill_octree->cubes_properties.size()) - 1);
+            num_lines += context.output_lines.size() + context.temp_lines.size();
         }
+        // Collect the lines.
+        std::vector<Line> lines;
+        lines.reserve(num_lines);
+        for (auto &context : contexts) {
+            append(lines, context.output_lines);
+            for (const Line &line : context.temp_lines)
+                if (line.a.x() != std::numeric_limits<coord_t>::max())
+                    lines.emplace_back(line);
+        }
+        // Convert lines to polylines.
+        all_polylines.reserve(lines.size());
+        std::transform(lines.begin(), lines.end(), std::back_inserter(all_polylines), [](const Line& l) { return Polyline{ l.a, l.b }; });
     }
 
+    // Crop all polylines
+    all_polylines = intersection_pl(std::move(all_polylines), to_polygons(expolygon));
+
+#ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
+    {
+        static int iRun = 0;
+        export_infill_lines_to_svg(expolygon, all_polylines, debug_out_path("FillAdaptive-initial-%d.svg", iRun++));
+    }
+#endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
+
     if (params.dont_connect)
-    {
-        // Crop all polylines
-        polylines_out = intersection_pl(all_polylines, to_polygons(expolygon));
-    }
-    else
-    {
-        // Crop all polylines
-        all_polylines = intersection_pl(all_polylines, to_polygons(expolygon));
-
+        append(polylines_out, std::move(all_polylines));
+    else {
         Polylines boundary_polylines;
         Polylines non_boundary_polylines;
         for (const Polyline &polyline : all_polylines)
-        {
             // connect_infill required all polylines to touch the boundary.
-            if(polyline.lines().size() == 1 && expolygon.has_boundary_point(polyline.lines().front().a) && expolygon.has_boundary_point(polyline.lines().front().b))
-            {
+            if (polyline.lines().size() == 1 && expolygon.has_boundary_point(polyline.lines().front().a) && expolygon.has_boundary_point(polyline.lines().front().b))
                 boundary_polylines.push_back(polyline);
-            }
             else
-            {
                 non_boundary_polylines.push_back(polyline);
-            }
-        }
 
-        if(!boundary_polylines.empty())
-        {
+        if (!boundary_polylines.empty()) {
             boundary_polylines = chain_polylines(boundary_polylines);
             FillAdaptive::connect_infill(std::move(boundary_polylines), expolygon, polylines_out, this->spacing, params);
         }
 
-        polylines_out.insert(polylines_out.end(), non_boundary_polylines.begin(), non_boundary_polylines.end());
+        append(polylines_out, std::move(non_boundary_polylines));
     }
 
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
+#ifdef ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT
     {
-        static int iRuna = 0;
-        BoundingBox bbox_svg = this->bounding_box;
-        {
-            ::Slic3r::SVG svg(debug_out_path("FillAdaptive-%d.svg", iRuna), bbox_svg);
-            for (const Polyline &polyline : polylines_out)
-            {
-                for (const Line &line : polyline.lines())
-                {
-                    Point from = line.a;
-                    Point to = line.b;
-                    Point diff = to - from;
-
-                    float shrink_length = scale_(0.4);
-                    float line_slope = (float)diff.y() / diff.x();
-                    float shrink_x = shrink_length / (float)std::sqrt(1.0 + (line_slope * line_slope));
-                    float shrink_y = line_slope * shrink_x;
-
-                    to.x() -= shrink_x;
-                    to.y() -= shrink_y;
-                    from.x() += shrink_x;
-                    from.y() += shrink_y;
-
-                    svg.draw(Line(from, to));
-                }
-            }
-        }
-
-        iRuna++;
+        static int iRun = 0;
+        export_infill_lines_to_svg(expolygon, polylines_out, debug_out_path("FillAdaptive-final-%d.svg", iRun ++));
     }
-#endif /* SLIC3R_DEBUG */
+#endif /* ADAPTIVE_CUBIC_INFILL_DEBUG_OUTPUT */
 }
 
-void FillAdaptive::generate_infill_lines(
-        FillAdaptive_Internal::Cube *cube,
-        double z_position,
-        const Vec3d &origin,
-        const Transform3d &rotation_matrix,
-        std::vector<Lines> &dir_lines_out,
-        const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties,
-        int depth)
+static double bbox_max_radius(const BoundingBoxf3 &bbox, const Vec3d &center)
 {
-    using namespace FillAdaptive_Internal;
-
-    if(cube == nullptr)
-    {
-        return;
-    }
-
-    Vec3d cube_center_tranformed = rotation_matrix * cube->center;
-    double z_diff = std::abs(z_position - cube_center_tranformed.z());
-
-    if (z_diff > cubes_properties[depth].height / 2)
-    {
-        return;
-    }
-
-    if (z_diff < cubes_properties[depth].line_z_distance)
-    {
-        Point from(
-                scale_((cubes_properties[depth].diagonal_length / 2) * (cubes_properties[depth].line_z_distance - z_diff) / cubes_properties[depth].line_z_distance),
-                scale_(cubes_properties[depth].line_xy_distance - ((z_position - (cube_center_tranformed.z() - cubes_properties[depth].line_z_distance)) / sqrt(2))));
-        Point to(-from.x(), from.y());
-        // Relative to cube center
-
-        double rotation_angle = (2.0 * M_PI) / 3.0;
-        for (Lines &lines : dir_lines_out)
-        {
-            Vec3d offset = cube_center_tranformed - (rotation_matrix * origin);
-            Point from_abs(from), to_abs(to);
-
-            from_abs.x() += int(scale_(offset.x()));
-            from_abs.y() += int(scale_(offset.y()));
-            to_abs.x() += int(scale_(offset.x()));
-            to_abs.y() += int(scale_(offset.y()));
-
-//            lines.emplace_back(from_abs, to_abs);
-            this->connect_lines(lines, Line(from_abs, to_abs));
-
-            from.rotate(rotation_angle);
-            to.rotate(rotation_angle);
-        }
-    }
-
-    for(const std::unique_ptr<Cube> &child : cube->children)
-    {
-        if(child != nullptr)
-        {
-            generate_infill_lines(child.get(), z_position, origin, rotation_matrix, dir_lines_out, cubes_properties, depth - 1);
-        }
-    }
+    const auto p = (bbox.min - center);
+    const auto s = bbox.size();
+    double r2max = 0.;
+    for (int i = 0; i < 8; ++ i)
+        r2max = std::max(r2max, (p + Vec3d(s.x() * double(i & 1), s.y() * double(i & 2), s.z() * double(i & 4))).squaredNorm());
+    return sqrt(r2max);
 }
 
-void FillAdaptive::connect_lines(Lines &lines, Line new_line)
+static std::vector<FillAdaptive_Internal::CubeProperties> make_cubes_properties(double max_cube_edge_length, double line_spacing)
 {
-    auto eps = int(scale_(0.10));
-    for (size_t i = 0; i < lines.size(); ++i)
+    max_cube_edge_length += EPSILON;
+
+    std::vector<FillAdaptive_Internal::CubeProperties> cubes_properties;
+    for (double edge_length = line_spacing * 2.;; edge_length *= 2.)
     {
-        if (std::abs(new_line.a.x() - lines[i].b.x()) < eps && std::abs(new_line.a.y() - lines[i].b.y()) < eps)
-        {
-            new_line.a = lines[i].a;
-            lines.erase(lines.begin() + i);
-            --i;
-            continue;
-        }
-
-        if (std::abs(new_line.b.x() - lines[i].a.x()) < eps && std::abs(new_line.b.y() - lines[i].a.y()) < eps)
-        {
-            new_line.b = lines[i].b;
-            lines.erase(lines.begin() + i);
-            --i;
-            continue;
-        }
-    }
-
-    lines.emplace_back(new_line.a, new_line.b);
-}
-
-std::unique_ptr<FillAdaptive_Internal::Octree> FillAdaptive::build_octree(
-    TriangleMesh &triangle_mesh,
-    coordf_t line_spacing,
-    const Vec3d &cube_center)
-{
-    using namespace FillAdaptive_Internal;
-
-    if(line_spacing <= 0 || std::isnan(line_spacing))
-    {
-        return nullptr;
-    }
-
-    Vec3d bb_size = triangle_mesh.bounding_box().size();
-    // The furthest point from the center of the bottom of the mesh bounding box.
-    double furthest_point = std::sqrt(((bb_size.x() * bb_size.x()) / 4.0) +
-                                      ((bb_size.y() * bb_size.y()) / 4.0) +
-                                      (bb_size.z() * bb_size.z()));
-    double max_cube_edge_length = furthest_point * 2;
-
-    std::vector<CubeProperties> cubes_properties;
-    for (double edge_length = (line_spacing * 2); edge_length < (max_cube_edge_length * 2); edge_length *= 2)
-    {
-        CubeProperties props{};
+        FillAdaptive_Internal::CubeProperties props{};
         props.edge_length = edge_length;
         props.height = edge_length * sqrt(3);
         props.diagonal_length = edge_length * sqrt(2);
         props.line_z_distance = edge_length / sqrt(3);
         props.line_xy_distance = edge_length / sqrt(6);
-        cubes_properties.push_back(props);
+        cubes_properties.emplace_back(props);
+        if (edge_length > max_cube_edge_length)
+            break;
     }
+    return cubes_properties;
+}
 
-    if (triangle_mesh.its.vertices.empty())
-    {
-        triangle_mesh.require_shared_vertices();
+static inline bool is_overhang_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &up)
+{
+    // Calculate triangle normal.
+    auto n = (b - a).cross(c - b);
+    return n.dot(up) > 0.707 * n.norm();
+}
+
+static void transform_center(FillAdaptive_Internal::Cube *current_cube, const Eigen::Matrix3d &rot)
+{
+#ifndef NDEBUG
+    current_cube->center_octree = current_cube->center;
+#endif // NDEBUG
+    current_cube->center = rot * current_cube->center;
+    for (auto *child : current_cube->children)
+        if (child)
+            transform_center(child, rot);
+}
+
+FillAdaptive_Internal::OctreePtr FillAdaptive_Internal::build_octree(const indexed_triangle_set &triangle_mesh, const Vec3d &up_vector, coordf_t line_spacing, bool support_overhangs_only)
+{
+    assert(line_spacing > 0);
+    assert(! std::isnan(line_spacing));
+
+    BoundingBox3Base<Vec3f>     bbox(triangle_mesh.vertices);
+    Vec3d                       cube_center      = bbox.center().cast<double>();
+    std::vector<CubeProperties> cubes_properties = make_cubes_properties(double(bbox.size().maxCoeff()), line_spacing);
+    auto                        octree           = OctreePtr(new Octree(cube_center, cubes_properties));
+
+    if (cubes_properties.size() > 1) {
+        for (auto &tri : triangle_mesh.indices) {
+            auto a = triangle_mesh.vertices[tri[0]].cast<double>();
+            auto b = triangle_mesh.vertices[tri[1]].cast<double>();
+            auto c = triangle_mesh.vertices[tri[2]].cast<double>();
+            if (support_overhangs_only && ! is_overhang_triangle(a, b, c, up_vector))
+                continue;
+            double edge_length_half = 0.5 * cubes_properties.back().edge_length;
+            Vec3d  diag_half(edge_length_half, edge_length_half, edge_length_half);
+            octree->insert_triangle(
+                a, b, c,
+                octree->root_cube, 
+                BoundingBoxf3(octree->root_cube->center - diag_half, octree->root_cube->center + diag_half),
+                int(cubes_properties.size()) - 1);
+        }
+        {
+            // Transform the octree to world coordinates to reduce computation when extracting infill lines.
+            auto rot = adaptive_fill_octree_transform_to_world().toRotationMatrix();
+            transform_center(octree->root_cube, rot);
+            octree->origin = rot * octree->origin;
+        }
     }
 
-    AABBTreeIndirect::Tree3f aabbTree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
-            triangle_mesh.its.vertices, triangle_mesh.its.indices);
-    auto octree = std::make_unique<Octree>(std::make_unique<Cube>(cube_center), cube_center, cubes_properties);
-
-    FillAdaptive::expand_cube(octree->root_cube.get(), cubes_properties, aabbTree, triangle_mesh, int(cubes_properties.size()) - 1);
-
     return octree;
 }
 
-void FillAdaptive::expand_cube(
-    FillAdaptive_Internal::Cube *cube,
-    const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties,
-    const AABBTreeIndirect::Tree3f &distance_tree,
-    const TriangleMesh &triangle_mesh, int depth)
+void FillAdaptive_Internal::Octree::insert_triangle(const Vec3d &a, const Vec3d &b, const Vec3d &c, Cube *current_cube, const BoundingBoxf3 &current_bbox, int depth)
 {
-    using namespace FillAdaptive_Internal;
+    assert(current_cube);
+    assert(depth > 0);
 
-    if (cube == nullptr || depth == 0)
-    {
-        return;
-    }
-
-    std::vector<Vec3d> child_centers = {
-        Vec3d(-1, -1, -1), Vec3d( 1, -1, -1), Vec3d(-1,  1, -1), Vec3d( 1,  1, -1),
-        Vec3d(-1, -1,  1), Vec3d( 1, -1,  1), Vec3d(-1,  1,  1), Vec3d( 1,  1,  1)
-    };
-
-    double cube_radius_squared = (cubes_properties[depth].height * cubes_properties[depth].height) / 16;
-
-    for (size_t i = 0; i < 8; ++i)
-    {
+    for (size_t i = 0; i < 8; ++ i) {
         const Vec3d &child_center = child_centers[i];
-        Vec3d child_center_transformed = cube->center + (child_center * (cubes_properties[depth].edge_length / 4));
-
-        if(AABBTreeIndirect::is_any_triangle_in_radius(triangle_mesh.its.vertices, triangle_mesh.its.indices,
-            distance_tree, child_center_transformed, cube_radius_squared))
-        {
-            cube->children[i] = std::make_unique<Cube>(child_center_transformed);
-            FillAdaptive::expand_cube(cube->children[i].get(), cubes_properties, distance_tree, triangle_mesh, depth - 1);
-        }
-    }
-}
-
-void FillAdaptive_Internal::Octree::propagate_point(
-    Vec3d                                                     point,
-    FillAdaptive_Internal::Cube *                             current,
-    int                                                       depth,
-    const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties)
-{
-    using namespace FillAdaptive_Internal;
-
-    if(depth <= 0)
-    {
-        return;
-    }
-
-    size_t octant_idx = Octree::find_octant(point, current->center);
-    Cube * child = current->children[octant_idx].get();
-
-    // Octant not exists, then create it
-    if(child == nullptr) {
-        std::vector<Vec3d> child_centers = {
-            Vec3d(-1, -1, -1), Vec3d( 1, -1, -1), Vec3d(-1,  1, -1), Vec3d( 1,  1, -1),
-            Vec3d(-1, -1,  1), Vec3d( 1, -1,  1), Vec3d(-1,  1,  1), Vec3d( 1,  1,  1)
-        };
-
-        const Vec3d &child_center = child_centers[octant_idx];
-        Vec3d child_center_transformed = current->center + (child_center * (cubes_properties[depth].edge_length / 4));
-
-        current->children[octant_idx] = std::make_unique<Cube>(child_center_transformed);
-        child = current->children[octant_idx].get();
-    }
-
-    Octree::propagate_point(point, child, (depth - 1), cubes_properties);
-}
-
-std::unique_ptr<FillAdaptive_Internal::Octree> FillSupportCubic::build_octree(
-    TriangleMesh &     triangle_mesh,
-    coordf_t           line_spacing,
-    const Vec3d &      cube_center,
-    const Transform3d &rotation_matrix)
-{
-    using namespace FillAdaptive_Internal;
-
-    if(line_spacing <= 0 || std::isnan(line_spacing))
-    {
-        return nullptr;
-    }
-
-    Vec3d bb_size = triangle_mesh.bounding_box().size();
-    // The furthest point from the center of the bottom of the mesh bounding box.
-    double furthest_point = std::sqrt(((bb_size.x() * bb_size.x()) / 4.0) +
-                                      ((bb_size.y() * bb_size.y()) / 4.0) +
-                                      (bb_size.z() * bb_size.z()));
-    double max_cube_edge_length = furthest_point * 2;
-
-    std::vector<CubeProperties> cubes_properties;
-    for (double edge_length = (line_spacing * 2); edge_length < (max_cube_edge_length * 2); edge_length *= 2)
-    {
-        CubeProperties props{};
-        props.edge_length = edge_length;
-        props.height = edge_length * sqrt(3);
-        props.diagonal_length = edge_length * sqrt(2);
-        props.line_z_distance = edge_length / sqrt(3);
-        props.line_xy_distance = edge_length / sqrt(6);
-        cubes_properties.push_back(props);
-    }
-
-    if (triangle_mesh.its.vertices.empty())
-    {
-        triangle_mesh.require_shared_vertices();
-    }
-
-    AABBTreeIndirect::Tree3f aabbTree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
-        triangle_mesh.its.vertices, triangle_mesh.its.indices);
-
-    auto octree = std::make_unique<Octree>(std::make_unique<Cube>(cube_center), cube_center, cubes_properties);
-
-    double cube_edge_length = line_spacing / 2.0;
-    int max_depth = int(octree->cubes_properties.size()) - 1;
-    BoundingBoxf3 mesh_bb = triangle_mesh.bounding_box();
-    Vec3f vertical(0, 0, 1);
-
-    for (size_t facet_idx = 0; facet_idx < triangle_mesh.stl.facet_start.size(); ++facet_idx)
-    {
-        if(triangle_mesh.stl.facet_start[facet_idx].normal.dot(vertical) <= 0.707)
-        {
-            // The angle is smaller than PI/4, than infill don't to be there
-            continue;
-        }
-
-        stl_vertex v_1 = triangle_mesh.stl.facet_start[facet_idx].vertex[0];
-        stl_vertex v_2 = triangle_mesh.stl.facet_start[facet_idx].vertex[1];
-        stl_vertex v_3 = triangle_mesh.stl.facet_start[facet_idx].vertex[2];
-
-        std::vector<Vec3d> triangle_vertices =
-            {Vec3d(v_1.x(), v_1.y(), v_1.z()),
-             Vec3d(v_2.x(), v_2.y(), v_2.z()),
-             Vec3d(v_3.x(), v_3.y(), v_3.z())};
-
-        BoundingBoxf3 triangle_bb(triangle_vertices);
-
-        Vec3d triangle_start_relative = triangle_bb.min - mesh_bb.min;
-        Vec3d triangle_end_relative   = triangle_bb.max - mesh_bb.min;
-
-        Vec3crd triangle_start_idx = Vec3crd(
-            int(std::floor(triangle_start_relative.x() / cube_edge_length)),
-            int(std::floor(triangle_start_relative.y() / cube_edge_length)),
-            int(std::floor(triangle_start_relative.z() / cube_edge_length)));
-        Vec3crd triangle_end_idx = Vec3crd(
-            int(std::floor(triangle_end_relative.x() / cube_edge_length)),
-            int(std::floor(triangle_end_relative.y() / cube_edge_length)),
-            int(std::floor(triangle_end_relative.z() / cube_edge_length)));
-
-        for (int z = triangle_start_idx.z(); z <= triangle_end_idx.z(); ++z)
-        {
-            for (int y = triangle_start_idx.y(); y <= triangle_end_idx.y(); ++y)
-            {
-                for (int x = triangle_start_idx.x(); x <= triangle_end_idx.x(); ++x)
-                {
-                    Vec3d cube_center_relative(x * cube_edge_length + (cube_edge_length / 2.0), y * cube_edge_length + (cube_edge_length / 2.0), z * cube_edge_length);
-                    Vec3d cube_center_absolute = cube_center_relative + mesh_bb.min;
-
-                    double cube_center_absolute_arr[3] = {cube_center_absolute.x(), cube_center_absolute.y(), cube_center_absolute.z()};
-                    double distance = 0, cord_u = 0, cord_v = 0;
-
-                    double dir[3] = {0.0, 0.0, 1.0};
-
-                    double vert_0[3] = {triangle_vertices[0].x(),
-                                        triangle_vertices[0].y(),
-                                        triangle_vertices[0].z()};
-                    double vert_1[3] = {triangle_vertices[1].x(),
-                                        triangle_vertices[1].y(),
-                                        triangle_vertices[1].z()};
-                    double vert_2[3] = {triangle_vertices[2].x(),
-                                        triangle_vertices[2].y(),
-                                        triangle_vertices[2].z()};
-
-                    if(intersect_triangle(cube_center_absolute_arr, dir, vert_0, vert_1, vert_2, &distance, &cord_u, &cord_v) && distance > 0 && distance <= cube_edge_length)
-                    {
-                        Vec3d cube_center_transformed(cube_center_absolute.x(), cube_center_absolute.y(), cube_center_absolute.z() + (cube_edge_length / 2.0));
-                        Octree::propagate_point(rotation_matrix * cube_center_transformed, octree->root_cube.get(), max_depth, octree->cubes_properties);
-                    }
-                }
+        // Calculate a slightly expanded bounding box of a child cube to cope with triangles touching a cube wall and other numeric errors.
+        // We will rather densify the octree a bit more than necessary instead of missing a triangle.
+        BoundingBoxf3 bbox;
+        for (int k = 0; k < 3; ++ k) {
+            if (child_center[k] == -1.) {
+                bbox.min[k] = current_bbox.min[k];
+                bbox.max[k] = current_cube->center[k] + EPSILON;
+            } else {
+                bbox.min[k] = current_cube->center[k] - EPSILON;
+                bbox.max[k] = current_bbox.max[k];
             }
         }
+        if (triangle_AABB_intersects(a, b, c, bbox)) {
+            if (! current_cube->children[i])
+                current_cube->children[i] = this->pool.construct(current_cube->center + (child_center * (this->cubes_properties[depth].edge_length / 4)));
+            if (depth > 1)
+                this->insert_triangle(a, b, c, current_cube->children[i], bbox, depth - 1);
+        }
     }
-
-    return octree;
-}
-
-void FillSupportCubic::_fill_surface_single(const FillParams &             params,
-                                            unsigned int                   thickness_layers,
-                                            const std::pair<float, Point> &direction,
-                                            ExPolygon &                    expolygon,
-                                            Polylines &                    polylines_out)
-{
-    if (this->support_fill_octree != nullptr)
-        this->generate_infill(params, thickness_layers, direction, expolygon, polylines_out, this->support_fill_octree);
 }
 
 } // namespace Slic3r
diff --git a/src/libslic3r/Fill/FillAdaptive.hpp b/src/libslic3r/Fill/FillAdaptive.hpp
index b24f206da..906414747 100644
--- a/src/libslic3r/Fill/FillAdaptive.hpp
+++ b/src/libslic3r/Fill/FillAdaptive.hpp
@@ -1,3 +1,13 @@
+// Adaptive cubic infill was inspired by the work of @mboerwinkle
+// as implemented for Cura.
+// https://github.com/Ultimaker/CuraEngine/issues/381
+// https://github.com/Ultimaker/CuraEngine/pull/401
+//
+// Our implementation is more accurate (discretizes a bit less cubes than Cura's)
+// by splitting only such cubes which contain a triangle. 
+// Our line extraction is time optimal instead of O(n^2) when connecting extracted lines,
+// and we also implemented adaptivity for supporting internal overhangs only.
+
 #ifndef slic3r_FillAdaptive_hpp_
 #define slic3r_FillAdaptive_hpp_
 
@@ -5,49 +15,39 @@
 
 #include "FillBase.hpp"
 
+struct indexed_triangle_set;
+
 namespace Slic3r {
 
 class PrintObject;
-class TriangleMesh;
-
 namespace FillAdaptive_Internal
 {
-    struct CubeProperties
-    {
-        double edge_length;     // Lenght of edge of a cube
-        double height;          // Height of rotated cube (standing on the corner)
-        double diagonal_length; // Length of diagonal of a cube a face
-        double line_z_distance; // Defines maximal distance from a center of a cube on Z axis on which lines will be created
-        double line_xy_distance;// Defines maximal distance from a center of a cube on X and Y axis on which lines will be created
+    struct Octree;
+    // To keep the definition of Octree opaque, we have to define a custom deleter.
+    struct OctreeDeleter {
+        void operator()(Octree *p);
     };
+    using OctreePtr = std::unique_ptr<Octree, OctreeDeleter>;
 
-    struct Cube
-    {
-        Vec3d center;
-        std::unique_ptr<Cube> children[8] = {};
-        Cube(const Vec3d &center) : center(center) {}
-    };
+    // Calculate line spacing for
+    // 1) adaptive cubic infill
+    // 2) adaptive internal support cubic infill
+    // Returns zero for a particular infill type if no such infill is to be generated.
+    std::pair<double, double>        adaptive_fill_line_spacing(const PrintObject &print_object);
 
-    struct Octree
-    {
-        std::unique_ptr<Cube> root_cube;
-        Vec3d origin;
-        std::vector<CubeProperties> cubes_properties;
+    // Rotation of the octree to stand on one of its corners.
+    Eigen::Quaterniond               adaptive_fill_octree_transform_to_world();
+    // Inverse roation of the above.
+    Eigen::Quaterniond               adaptive_fill_octree_transform_to_octree();
 
-        Octree(std::unique_ptr<Cube> rootCube, const Vec3d &origin, const std::vector<CubeProperties> &cubes_properties)
-            : root_cube(std::move(rootCube)), origin(origin), cubes_properties(cubes_properties) {}
-
-        inline static int find_octant(const Vec3d &i_cube, const Vec3d &current)
-        {
-            return (i_cube.z() > current.z()) * 4 + (i_cube.y() > current.y()) * 2 + (i_cube.x() > current.x());
-        }
-
-        static void propagate_point(
-            Vec3d                        point,
-            FillAdaptive_Internal::Cube *current_cube,
-            int                          depth,
-            const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties);
-    };
+    FillAdaptive_Internal::OctreePtr build_octree(
+        // Mesh is rotated to the coordinate system of the octree.
+        const indexed_triangle_set  &triangle_mesh, 
+        // Up vector of the mesh rotated to the coordinate system of the octree.
+        const Vec3d                 &up_vector, 
+        coordf_t                     line_spacing, 
+        // If true, octree is densified below internal overhangs only.
+        bool                         support_overhangs_only);
 }; // namespace FillAdaptive_Internal
 
 //
@@ -70,70 +70,8 @@ protected:
 	    Polylines                       &polylines_out);
 
 	virtual bool no_sort() const { return true; }
-
-    void generate_infill_lines(
-        FillAdaptive_Internal::Cube *cube,
-        double                       z_position,
-        const Vec3d &                origin,
-        const Transform3d &          rotation_matrix,
-        std::vector<Lines> &         dir_lines_out,
-        const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties,
-        int  depth);
-
-    static void connect_lines(Lines &lines, Line new_line);
-
-    void generate_infill(const FillParams &             params,
-                         unsigned int                   thickness_layers,
-                         const std::pair<float, Point> &direction,
-                         ExPolygon &                    expolygon,
-                         Polylines &                    polylines_out,
-                         FillAdaptive_Internal::Octree *octree);
-
-public:
-    static std::unique_ptr<FillAdaptive_Internal::Octree> build_octree(
-        TriangleMesh &triangle_mesh,
-        coordf_t      line_spacing,
-        const Vec3d & cube_center);
-
-    static void expand_cube(
-        FillAdaptive_Internal::Cube *cube,
-        const std::vector<FillAdaptive_Internal::CubeProperties> &cubes_properties,
-        const AABBTreeIndirect::Tree3f &distance_tree,
-        const TriangleMesh &            triangle_mesh,
-        int                             depth);
 };
 
-class FillSupportCubic : public FillAdaptive
-{
-public:
-    virtual ~FillSupportCubic() = default;
-
-protected:
-    virtual Fill* clone() const { return new FillSupportCubic(*this); };
-
-    virtual bool no_sort() const { return true; }
-
-    virtual void _fill_surface_single(
-        const FillParams                &params,
-        unsigned int                     thickness_layers,
-        const std::pair<float, Point>   &direction,
-        ExPolygon                       &expolygon,
-        Polylines                       &polylines_out);
-
-public:
-    static std::unique_ptr<FillAdaptive_Internal::Octree> build_octree(
-        TriangleMesh &     triangle_mesh,
-        coordf_t           line_spacing,
-        const Vec3d &      cube_center,
-        const Transform3d &rotation_matrix);
-};
-
-// Calculate line spacing for
-// 1) adaptive cubic infill
-// 2) adaptive internal support cubic infill
-// Returns zero for a particular infill type if no such infill is to be generated.
-std::pair<double, double> adaptive_fill_line_spacing(const PrintObject &print_object);
-
 } // namespace Slic3r
 
 #endif // slic3r_FillAdaptive_hpp_
diff --git a/src/libslic3r/Fill/FillBase.cpp b/src/libslic3r/Fill/FillBase.cpp
index 077555d2c..5866330b9 100644
--- a/src/libslic3r/Fill/FillBase.cpp
+++ b/src/libslic3r/Fill/FillBase.cpp
@@ -39,7 +39,7 @@ Fill* Fill::new_from_type(const InfillPattern type)
     case ipHilbertCurve:        return new FillHilbertCurve();
     case ipOctagramSpiral:      return new FillOctagramSpiral();
     case ipAdaptiveCubic:       return new FillAdaptive();
-    case ipSupportCubic:        return new FillSupportCubic();
+    case ipSupportCubic:        return new FillAdaptive();
     default: throw Slic3r::InvalidArgument("unknown type");
     }
 }
diff --git a/src/libslic3r/Fill/FillBase.hpp b/src/libslic3r/Fill/FillBase.hpp
index 77620e118..4d822ddea 100644
--- a/src/libslic3r/Fill/FillBase.hpp
+++ b/src/libslic3r/Fill/FillBase.hpp
@@ -77,8 +77,6 @@ public:
 
     // Octree builds on mesh for usage in the adaptive cubic infill
     FillAdaptive_Internal::Octree* adapt_fill_octree = nullptr;
-    // Octree builds on mesh for usage in the support cubic infill
-    FillAdaptive_Internal::Octree* support_fill_octree = nullptr;
 
 public:
     virtual ~Fill() {}
diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp
index 75f3708d2..b690b478d 100644
--- a/src/libslic3r/Geometry.hpp
+++ b/src/libslic3r/Geometry.hpp
@@ -281,7 +281,7 @@ bool directions_parallel(double angle1, double angle2, double max_diff = 0);
 template<class T> bool contains(const std::vector<T> &vector, const Point &point);
 template<typename T> T rad2deg(T angle) { return T(180.0) * angle / T(PI); }
 double rad2deg_dir(double angle);
-template<typename T> T deg2rad(T angle) { return T(PI) * angle / T(180.0); }
+template<typename T> constexpr T deg2rad(const T angle) { return T(PI) * angle / T(180.0); }
 template<typename T> T angle_to_0_2PI(T angle)
 {
     static const T TWO_PI = T(2) * T(PI);
diff --git a/src/libslic3r/Model.cpp b/src/libslic3r/Model.cpp
index dc06e3102..3539cc0fa 100644
--- a/src/libslic3r/Model.cpp
+++ b/src/libslic3r/Model.cpp
@@ -777,6 +777,38 @@ TriangleMesh ModelObject::raw_mesh() const
     return mesh;
 }
 
+// Non-transformed (non-rotated, non-scaled, non-translated) sum of non-modifier object volumes.
+// Currently used by ModelObject::mesh(), to calculate the 2D envelope for 2D plater
+// and to display the object statistics at ModelObject::print_info().
+indexed_triangle_set ModelObject::raw_indexed_triangle_set() const
+{
+    size_t num_vertices = 0;
+    size_t num_faces    = 0;
+    for (const ModelVolume *v : this->volumes)
+        if (v->is_model_part()) {
+            num_vertices += v->mesh().its.vertices.size();
+            num_faces    += v->mesh().its.indices.size();
+        }
+    indexed_triangle_set out;
+    out.vertices.reserve(num_vertices);
+    out.indices.reserve(num_faces);
+    for (const ModelVolume *v : this->volumes)
+        if (v->is_model_part()) {
+            size_t i = out.vertices.size();
+            size_t j = out.indices.size();
+            append(out.vertices, v->mesh().its.vertices);
+            append(out.indices,  v->mesh().its.indices);
+            auto m = v->get_matrix();
+            for (; i < out.vertices.size(); ++ i)
+                out.vertices[i] = (m * out.vertices[i].cast<double>()).cast<float>().eval();
+            if (v->is_left_handed()) {
+                for (; j < out.indices.size(); ++ j)
+                    std::swap(out.indices[j][0], out.indices[j][1]);
+            }
+        }
+    return out;
+}
+
 // Non-transformed (non-rotated, non-scaled, non-translated) sum of all object volumes.
 TriangleMesh ModelObject::full_raw_mesh() const
 {
diff --git a/src/libslic3r/Model.hpp b/src/libslic3r/Model.hpp
index a623f5cca..b9045e28b 100644
--- a/src/libslic3r/Model.hpp
+++ b/src/libslic3r/Model.hpp
@@ -244,6 +244,8 @@ public:
     // Non-transformed (non-rotated, non-scaled, non-translated) sum of non-modifier object volumes.
     // Currently used by ModelObject::mesh() and to calculate the 2D envelope for 2D plater.
     TriangleMesh raw_mesh() const;
+    // The same as above, but producing a lightweight indexed_triangle_set.
+    indexed_triangle_set raw_indexed_triangle_set() const;
     // Non-transformed (non-rotated, non-scaled, non-translated) sum of all object volumes.
     TriangleMesh full_raw_mesh() const;
     // A transformed snug bounding box around the non-modifier object volumes, without the translation applied.
diff --git a/src/libslic3r/Point.cpp b/src/libslic3r/Point.cpp
index c2417d0dc..f3ed41342 100644
--- a/src/libslic3r/Point.cpp
+++ b/src/libslic3r/Point.cpp
@@ -44,16 +44,6 @@ Pointf3s transform(const Pointf3s& points, const Transform3d& t)
     return ret_points;
 }
 
-void Point::rotate(double angle)
-{
-    double cur_x = (double)(*this)(0);
-    double cur_y = (double)(*this)(1);
-    double s     = ::sin(angle);
-    double c     = ::cos(angle);
-    (*this)(0) = (coord_t)round(c * cur_x - s * cur_y);
-    (*this)(1) = (coord_t)round(c * cur_y + s * cur_x);
-}
-
 void Point::rotate(double angle, const Point &center)
 {
     double cur_x = (double)(*this)(0);
diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp
index 30a1a4942..5082bb746 100644
--- a/src/libslic3r/Point.hpp
+++ b/src/libslic3r/Point.hpp
@@ -105,6 +105,7 @@ public:
     template<typename OtherDerived>
     Point(const Eigen::MatrixBase<OtherDerived> &other) : Vec2crd(other) {}
     static Point new_scale(coordf_t x, coordf_t y) { return Point(coord_t(scale_(x)), coord_t(scale_(y))); }
+    static Point new_scale(const Vec2d &v) { return Point(coord_t(scale_(v.x())), coord_t(scale_(v.y()))); }
 
     // This method allows you to assign Eigen expressions to MyVectorType
     template<typename OtherDerived>
@@ -121,7 +122,14 @@ public:
 	Point& operator*=(const double &rhs) { (*this)(0) = coord_t((*this)(0) * rhs); (*this)(1) = coord_t((*this)(1) * rhs); return *this; }
     Point operator*(const double &rhs) { return Point((*this)(0) * rhs, (*this)(1) * rhs); }
 
-    void   rotate(double angle);
+    void   rotate(double angle) { this->rotate(std::cos(angle), std::sin(angle)); }
+    void   rotate(double cos_a, double sin_a) {
+        double cur_x = (double)(*this)(0);
+        double cur_y = (double)(*this)(1);
+        (*this)(0) = (coord_t)round(cos_a * cur_x - sin_a * cur_y);
+        (*this)(1) = (coord_t)round(cos_a * cur_y + sin_a * cur_x);
+    }
+
     void   rotate(double angle, const Point &center);
     Point  rotated(double angle) const { Point res(*this); res.rotate(angle); return res; }
     Point  rotated(double angle, const Point &center) const { Point res(*this); res.rotate(angle, center); return res; }
diff --git a/src/libslic3r/Print.hpp b/src/libslic3r/Print.hpp
index 98a131411..471484005 100644
--- a/src/libslic3r/Print.hpp
+++ b/src/libslic3r/Print.hpp
@@ -32,6 +32,8 @@ class SupportLayer;
 
 namespace FillAdaptive_Internal {
     struct Octree;
+    struct OctreeDeleter;
+    using OctreePtr = std::unique_ptr<Octree, OctreeDeleter>;
 };
 
 // Print step IDs for keeping track of the print state.
@@ -239,7 +241,7 @@ private:
     void discover_horizontal_shells();
     void combine_infill();
     void _generate_support_material();
-    std::pair<std::unique_ptr<FillAdaptive_Internal::Octree>, std::unique_ptr<FillAdaptive_Internal::Octree>> prepare_adaptive_infill_data();
+    std::pair<FillAdaptive_Internal::OctreePtr, FillAdaptive_Internal::OctreePtr> prepare_adaptive_infill_data();
 
     // XYZ in scaled coordinates
     Vec3crd									m_size;
diff --git a/src/libslic3r/PrintObject.cpp b/src/libslic3r/PrintObject.cpp
index ddd41af01..0968f6cfc 100644
--- a/src/libslic3r/PrintObject.cpp
+++ b/src/libslic3r/PrintObject.cpp
@@ -434,74 +434,27 @@ void PrintObject::generate_support_material()
     }
 }
 
-//#define ADAPTIVE_SUPPORT_SIMPLE
-
-std::pair<std::unique_ptr<FillAdaptive_Internal::Octree>, std::unique_ptr<FillAdaptive_Internal::Octree>> PrintObject::prepare_adaptive_infill_data()
+std::pair<FillAdaptive_Internal::OctreePtr, FillAdaptive_Internal::OctreePtr> PrintObject::prepare_adaptive_infill_data()
 {
     using namespace FillAdaptive_Internal;
 
     auto [adaptive_line_spacing, support_line_spacing] = adaptive_fill_line_spacing(*this);
-
-    std::unique_ptr<Octree> adaptive_fill_octree = {}, support_fill_octree = {};
-
     if (adaptive_line_spacing == 0. && support_line_spacing == 0.)
-        return std::make_pair(std::move(adaptive_fill_octree), std::move(support_fill_octree));
+        return std::make_pair(OctreePtr(), OctreePtr());
 
-    TriangleMesh mesh = this->model_object()->raw_mesh();
-    mesh.transform(m_trafo, true);
-    // Apply XY shift
-    mesh.translate(- unscale<float>(m_center_offset.x()), - unscale<float>(m_center_offset.y()), 0);
-    // Center of the first cube in octree
-    Vec3d mesh_origin = mesh.bounding_box().center();
-
-#ifdef ADAPTIVE_SUPPORT_SIMPLE
-    if (mesh.its.vertices.empty())
+    indexed_triangle_set mesh = this->model_object()->raw_indexed_triangle_set();
+    Vec3d                up;
     {
-        mesh.require_shared_vertices();
-    }
-
-    Vec3f vertical(0, 0, 1);
-
-    indexed_triangle_set its_set;
-    its_set.vertices = mesh.its.vertices;
-
-    // Filter out non overhanging faces
-    for (size_t i = 0; i < mesh.its.indices.size(); ++i) {
-        stl_triangle_vertex_indices vertex_idx = mesh.its.indices[i];
-
-        auto its_calculate_normal = [](const stl_triangle_vertex_indices &index, const std::vector<stl_vertex> &vertices) {
-            stl_normal normal = (vertices[index.y()] - vertices[index.x()]).cross(vertices[index.z()] - vertices[index.x()]);
-            return normal;
-        };
-
-        stl_normal normal = its_calculate_normal(vertex_idx, mesh.its.vertices);
-        stl_normalize_vector(normal);
-
-        if(normal.dot(vertical) >= 0.707) {
-            its_set.indices.push_back(vertex_idx);
-        }
-    }
-
-    mesh = TriangleMesh(its_set);
-
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-    Slic3r::store_stl(debug_out_path("overhangs.stl").c_str(), &mesh, false);
-#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-#endif /* ADAPTIVE_SUPPORT_SIMPLE */
-
-    Vec3d rotation = Vec3d((5.0 * M_PI) / 4.0, Geometry::deg2rad(215.264), M_PI / 6.0);
-    Transform3d rotation_matrix = Geometry::assemble_transform(Vec3d::Zero(), rotation, Vec3d::Ones(), Vec3d::Ones()).inverse();
-
-    if (adaptive_line_spacing != 0.) {
+        auto m = adaptive_fill_octree_transform_to_octree().toRotationMatrix();
+        up = m * Vec3d(0., 0., 1.);
         // Rotate mesh and build octree on it with axis-aligned (standart base) cubes
-        mesh.transform(rotation_matrix);
-        adaptive_fill_octree = FillAdaptive::build_octree(mesh, adaptive_line_spacing, rotation_matrix * mesh_origin);
+        Transform3d m2 = m_trafo;
+        m2.translate(Vec3d(- unscale<float>(m_center_offset.x()), - unscale<float>(m_center_offset.y()), 0));
+        its_transform(mesh, m * m2, true);
     }
-
-    if (support_line_spacing != 0.)
-        support_fill_octree = FillSupportCubic::build_octree(mesh, support_line_spacing, rotation_matrix * mesh_origin, rotation_matrix);
-
-    return std::make_pair(std::move(adaptive_fill_octree), std::move(support_fill_octree));
+    return std::make_pair(
+        adaptive_line_spacing ? build_octree(mesh, up, adaptive_line_spacing, false) : OctreePtr(),
+        support_line_spacing  ? build_octree(mesh, up, support_line_spacing, true) : OctreePtr());
 }
 
 void PrintObject::clear_layers()