diff --git a/src/admesh/stl.h b/src/admesh/stl.h
index e0f2865f0..f7b8937ad 100644
--- a/src/admesh/stl.h
+++ b/src/admesh/stl.h
@@ -135,7 +135,7 @@ struct stl_file {
 	std::vector<stl_facet>     		facet_start;
 	std::vector<stl_neighbors> 		neighbors_start;
 	// Statistics
-	stl_stats     					stats;
+	stl_stats     					      stats;
 };
 
 struct indexed_triangle_set
@@ -149,9 +149,9 @@ struct indexed_triangle_set
 	}
 
 	std::vector<stl_triangle_vertex_indices> 	indices;
-	std::vector<stl_vertex>       				vertices;
+	std::vector<stl_vertex>       				    vertices;
 	//FIXME add normals once we get rid of the stl_file from TriangleMesh completely.
-	//std::vector<stl_normal> 					normals
+	//std::vector<stl_normal> 					      normals
 };
 
 extern bool stl_open(stl_file *stl, const char *file);
diff --git a/src/clipper/clipper.cpp b/src/clipper/clipper.cpp
index 0285d9167..e46a0797c 100644
--- a/src/clipper/clipper.cpp
+++ b/src/clipper/clipper.cpp
@@ -105,6 +105,15 @@ struct OutRec {
 
 //------------------------------------------------------------------------------
 
+inline IntPoint IntPoint2d(cInt x, cInt y)
+{
+  return IntPoint(x, y
+#ifdef CLIPPERLIB_USE_XYZ
+    , 0
+#endif // CLIPPERLIB_USE_XYZ
+    );
+}
+
 inline cInt Round(double val)
 {
   return static_cast<cInt>((val < 0) ? (val - 0.5) : (val + 0.5));
@@ -2243,7 +2252,7 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge)
                 while (maxIt != m_Maxima.end() && *maxIt < e->Curr.x()) 
                 {
                   if (horzEdge->OutIdx >= 0 && !IsOpen)
-                    AddOutPt(horzEdge, IntPoint(*maxIt, horzEdge->Bot.y()));
+                    AddOutPt(horzEdge, IntPoint2d(*maxIt, horzEdge->Bot.y()));
                   ++maxIt;
                 }
             }
@@ -2252,7 +2261,7 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge)
                 while (maxRit != m_Maxima.rend() && *maxRit > e->Curr.x())
                 {
                   if (horzEdge->OutIdx >= 0 && !IsOpen)
-                    AddOutPt(horzEdge, IntPoint(*maxRit, horzEdge->Bot.y()));
+                    AddOutPt(horzEdge, IntPoint2d(*maxRit, horzEdge->Bot.y()));
                   ++maxRit;
                 }
             }
@@ -2297,12 +2306,12 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge)
         
 		if(dir == dLeftToRight)
         {
-          IntPoint Pt = IntPoint(e->Curr.x(), horzEdge->Curr.y());
+          IntPoint Pt = IntPoint2d(e->Curr.x(), horzEdge->Curr.y());
           IntersectEdges(horzEdge, e, Pt);
         }
         else
         {
-          IntPoint Pt = IntPoint(e->Curr.x(), horzEdge->Curr.y());
+          IntPoint Pt = IntPoint2d(e->Curr.x(), horzEdge->Curr.y());
           IntersectEdges( e, horzEdge, Pt);
         }
         TEdge* eNext = (dir == dLeftToRight) ? e->NextInAEL : e->PrevInAEL;
@@ -3372,14 +3381,14 @@ void ClipperOffset::AddPath(const Path& path, JoinType joinType, EndType endType
   //if this path's lowest pt is lower than all the others then update m_lowest
   if (endType != etClosedPolygon) return;
   if (m_lowest.x() < 0)
-    m_lowest = IntPoint(m_polyNodes.ChildCount() - 1, k);
+    m_lowest = IntPoint2d(m_polyNodes.ChildCount() - 1, k);
   else
   {
     IntPoint ip = m_polyNodes.Childs[(int)m_lowest.x()]->Contour[(int)m_lowest.y()];
     if (newNode->Contour[k].y() > ip.y() ||
       (newNode->Contour[k].y() == ip.y() &&
       newNode->Contour[k].x() < ip.x()))
-      m_lowest = IntPoint(m_polyNodes.ChildCount() - 1, k);
+      m_lowest = IntPoint2d(m_polyNodes.ChildCount() - 1, k);
   }
 }
 //------------------------------------------------------------------------------
@@ -3427,10 +3436,10 @@ void ClipperOffset::Execute(Paths& solution, double delta)
   {
     IntRect r = clpr.GetBounds();
     Path outer(4);
-    outer[0] = IntPoint(r.left - 10, r.bottom + 10);
-    outer[1] = IntPoint(r.right + 10, r.bottom + 10);
-    outer[2] = IntPoint(r.right + 10, r.top - 10);
-    outer[3] = IntPoint(r.left - 10, r.top - 10);
+    outer[0] = IntPoint2d(r.left - 10, r.bottom + 10);
+    outer[1] = IntPoint2d(r.right + 10, r.bottom + 10);
+    outer[2] = IntPoint2d(r.right + 10, r.top - 10);
+    outer[3] = IntPoint2d(r.left - 10, r.top - 10);
 
     clpr.AddPath(outer, ptSubject, true);
     clpr.ReverseSolution(true);
@@ -3457,10 +3466,10 @@ void ClipperOffset::Execute(PolyTree& solution, double delta)
   {
     IntRect r = clpr.GetBounds();
     Path outer(4);
-    outer[0] = IntPoint(r.left - 10, r.bottom + 10);
-    outer[1] = IntPoint(r.right + 10, r.bottom + 10);
-    outer[2] = IntPoint(r.right + 10, r.top - 10);
-    outer[3] = IntPoint(r.left - 10, r.top - 10);
+    outer[0] = IntPoint2d(r.left - 10, r.bottom + 10);
+    outer[1] = IntPoint2d(r.right + 10, r.bottom + 10);
+    outer[2] = IntPoint2d(r.right + 10, r.top - 10);
+    outer[3] = IntPoint2d(r.left - 10, r.top - 10);
 
     clpr.AddPath(outer, ptSubject, true);
     clpr.ReverseSolution(true);
@@ -3536,7 +3545,7 @@ void ClipperOffset::DoOffset(double delta)
         double X = 1.0, Y = 0.0;
         for (cInt j = 1; j <= steps; j++)
         {
-          m_destPoly.push_back(IntPoint(
+          m_destPoly.push_back(IntPoint2d(
             Round(m_srcPoly[0].x() + X * delta),
             Round(m_srcPoly[0].y() + Y * delta)));
           double X2 = X;
@@ -3549,7 +3558,7 @@ void ClipperOffset::DoOffset(double delta)
         double X = -1.0, Y = -1.0;
         for (int j = 0; j < 4; ++j)
         {
-          m_destPoly.push_back(IntPoint(
+          m_destPoly.push_back(IntPoint2d(
             Round(m_srcPoly[0].x() + X * delta),
             Round(m_srcPoly[0].y() + Y * delta)));
           if (X < 0) X = 1;
@@ -3604,9 +3613,9 @@ void ClipperOffset::DoOffset(double delta)
       if (node.m_endtype == etOpenButt)
       {
         int j = len - 1;
-        pt1 = IntPoint(Round(m_srcPoly[j].x() + m_normals[j].x() * delta), Round(m_srcPoly[j].y() + m_normals[j].y() * delta));
+        pt1 = IntPoint2d(Round(m_srcPoly[j].x() + m_normals[j].x() * delta), Round(m_srcPoly[j].y() + m_normals[j].y() * delta));
         m_destPoly.push_back(pt1);
-        pt1 = IntPoint(Round(m_srcPoly[j].x() - m_normals[j].x() * delta), Round(m_srcPoly[j].y() - m_normals[j].y() * delta));
+        pt1 = IntPoint2d(Round(m_srcPoly[j].x() - m_normals[j].x() * delta), Round(m_srcPoly[j].y() - m_normals[j].y() * delta));
         m_destPoly.push_back(pt1);
       }
       else
@@ -3631,9 +3640,9 @@ void ClipperOffset::DoOffset(double delta)
 
       if (node.m_endtype == etOpenButt)
       {
-        pt1 = IntPoint(Round(m_srcPoly[0].x() - m_normals[0].x() * delta), Round(m_srcPoly[0].y() - m_normals[0].y() * delta));
+        pt1 = IntPoint2d(Round(m_srcPoly[0].x() - m_normals[0].x() * delta), Round(m_srcPoly[0].y() - m_normals[0].y() * delta));
         m_destPoly.push_back(pt1);
-        pt1 = IntPoint(Round(m_srcPoly[0].x() + m_normals[0].x() * delta), Round(m_srcPoly[0].y() + m_normals[0].y() * delta));
+        pt1 = IntPoint2d(Round(m_srcPoly[0].x() + m_normals[0].x() * delta), Round(m_srcPoly[0].y() + m_normals[0].y() * delta));
         m_destPoly.push_back(pt1);
       }
       else
@@ -3661,7 +3670,7 @@ void ClipperOffset::OffsetPoint(int j, int& k, JoinType jointype)
     double cosA = (m_normals[k].x() * m_normals[j].x() + m_normals[j].y() * m_normals[k].y() ); 
     if (cosA > 0) // angle => 0 degrees
     {
-      m_destPoly.push_back(IntPoint(Round(m_srcPoly[j].x() + m_normals[k].x() * m_delta),
+      m_destPoly.push_back(IntPoint2d(Round(m_srcPoly[j].x() + m_normals[k].x() * m_delta),
         Round(m_srcPoly[j].y() + m_normals[k].y() * m_delta)));
       return; 
     }
@@ -3672,10 +3681,10 @@ void ClipperOffset::OffsetPoint(int j, int& k, JoinType jointype)
 
   if (m_sinA * m_delta < 0)
   {
-    m_destPoly.push_back(IntPoint(Round(m_srcPoly[j].x() + m_normals[k].x() * m_delta),
+    m_destPoly.push_back(IntPoint2d(Round(m_srcPoly[j].x() + m_normals[k].x() * m_delta),
       Round(m_srcPoly[j].y() + m_normals[k].y() * m_delta)));
     m_destPoly.push_back(m_srcPoly[j]);
-    m_destPoly.push_back(IntPoint(Round(m_srcPoly[j].x() + m_normals[j].x() * m_delta),
+    m_destPoly.push_back(IntPoint2d(Round(m_srcPoly[j].x() + m_normals[j].x() * m_delta),
       Round(m_srcPoly[j].y() + m_normals[j].y() * m_delta)));
   }
   else
@@ -3699,10 +3708,10 @@ void ClipperOffset::DoSquare(int j, int k)
 {
   double dx = std::tan(std::atan2(m_sinA,
       m_normals[k].x() * m_normals[j].x() + m_normals[k].y() * m_normals[j].y()) / 4);
-  m_destPoly.push_back(IntPoint(
+  m_destPoly.push_back(IntPoint2d(
       Round(m_srcPoly[j].x() + m_delta * (m_normals[k].x() - m_normals[k].y() * dx)),
       Round(m_srcPoly[j].y() + m_delta * (m_normals[k].y() + m_normals[k].x() * dx))));
-  m_destPoly.push_back(IntPoint(
+  m_destPoly.push_back(IntPoint2d(
       Round(m_srcPoly[j].x() + m_delta * (m_normals[j].x() + m_normals[j].y() * dx)),
       Round(m_srcPoly[j].y() + m_delta * (m_normals[j].y() - m_normals[j].x() * dx))));
 }
@@ -3711,7 +3720,7 @@ void ClipperOffset::DoSquare(int j, int k)
 void ClipperOffset::DoMiter(int j, int k, double r)
 {
   double q = m_delta / r;
-  m_destPoly.push_back(IntPoint(Round(m_srcPoly[j].x() + (m_normals[k].x() + m_normals[j].x()) * q),
+  m_destPoly.push_back(IntPoint2d(Round(m_srcPoly[j].x() + (m_normals[k].x() + m_normals[j].x()) * q),
       Round(m_srcPoly[j].y() + (m_normals[k].y() + m_normals[j].y()) * q)));
 }
 //------------------------------------------------------------------------------
@@ -3725,14 +3734,14 @@ void ClipperOffset::DoRound(int j, int k)
   double X = m_normals[k].x(), Y = m_normals[k].y(), X2;
   for (int i = 0; i < steps; ++i)
   {
-    m_destPoly.push_back(IntPoint(
+    m_destPoly.push_back(IntPoint2d(
         Round(m_srcPoly[j].x() + X * m_delta),
         Round(m_srcPoly[j].y() + Y * m_delta)));
     X2 = X;
     X = X * m_cos - m_sin * Y;
     Y = X2 * m_sin + Y * m_cos;
   }
-  m_destPoly.push_back(IntPoint(
+  m_destPoly.push_back(IntPoint2d(
   Round(m_srcPoly[j].x() + m_normals[j].x() * m_delta),
   Round(m_srcPoly[j].y() + m_normals[j].y() * m_delta)));
 }
@@ -4001,7 +4010,7 @@ void Minkowski(const Path& poly, const Path& path,
       Path p;
       p.reserve(polyCnt);
       for (size_t j = 0; j < poly.size(); ++j)
-        p.push_back(IntPoint(path[i].x() + poly[j].x(), path[i].y() + poly[j].y()));
+        p.push_back(IntPoint2d(path[i].x() + poly[j].x(), path[i].y() + poly[j].y()));
       pp.push_back(p);
     }
   else
@@ -4010,7 +4019,7 @@ void Minkowski(const Path& poly, const Path& path,
       Path p;
       p.reserve(polyCnt);
       for (size_t j = 0; j < poly.size(); ++j)
-        p.push_back(IntPoint(path[i].x() - poly[j].x(), path[i].y() - poly[j].y()));
+        p.push_back(IntPoint2d(path[i].x() - poly[j].x(), path[i].y() - poly[j].y()));
       pp.push_back(p);
     }
 
@@ -4045,7 +4054,7 @@ void TranslatePath(const Path& input, Path& output, const IntPoint& delta)
   //precondition: input != output
   output.resize(input.size());
   for (size_t i = 0; i < input.size(); ++i)
-    output[i] = IntPoint(input[i].x() + delta.x(), input[i].y() + delta.y());
+    output[i] = IntPoint2d(input[i].x() + delta.x(), input[i].y() + delta.y());
 }
 //------------------------------------------------------------------------------
 
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 0cf67597a..d98d0e10f 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -198,6 +198,8 @@ add_library(libslic3r STATIC
     Tesselate.hpp
     TriangleMesh.cpp
     TriangleMesh.hpp
+    TriangleMeshSlicer.cpp
+    TriangleMeshSlicer.hpp
     TriangulateWall.hpp
     TriangulateWall.cpp
     utils.cpp
diff --git a/src/libslic3r/Model.cpp b/src/libslic3r/Model.cpp
index 8d73782a4..4ddcbac39 100644
--- a/src/libslic3r/Model.cpp
+++ b/src/libslic3r/Model.cpp
@@ -3,6 +3,7 @@
 #include "ModelArrange.hpp"
 #include "Geometry.hpp"
 #include "MTUtils.hpp"
+#include "TriangleMeshSlicer.hpp"
 #include "TriangleSelector.hpp"
 
 #include "Format/AMF.hpp"
diff --git a/src/libslic3r/PrintObject.cpp b/src/libslic3r/PrintObject.cpp
index 46069849b..b77b13345 100644
--- a/src/libslic3r/PrintObject.cpp
+++ b/src/libslic3r/PrintObject.cpp
@@ -10,6 +10,7 @@
 #include "Surface.hpp"
 #include "Slicing.hpp"
 #include "Tesselate.hpp"
+#include "TriangleMeshSlicer.hpp"
 #include "Utils.hpp"
 #include "Fill/FillAdaptive.hpp"
 #include "Format/STL.hpp"
@@ -2221,9 +2222,8 @@ std::vector<ExPolygons> PrintObject::slice_volumes(
             auto callback = TriangleMeshSlicer::throw_on_cancel_callback_type([print](){print->throw_if_canceled();});
             // TriangleMeshSlicer needs shared vertices, also this calls the repair() function.
             mesh.require_shared_vertices();
-            TriangleMeshSlicer mslicer;
-            mslicer.init(&mesh, callback);
-			mslicer.slice(z, mode, slicing_mode_normal_below_layer, mode_below, float(m_config.slice_closing_radius.value), &layers, callback);
+            MeshSlicingParamsExtended params { { mode, slicing_mode_normal_below_layer, mode_below }, float(m_config.slice_closing_radius.value) };
+			slice_mesh(mesh, z, params, layers, callback);
             m_print->throw_if_canceled();
         }
     }
@@ -2245,13 +2245,14 @@ std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, S
 	        // apply XY shift
 	        mesh.translate(- unscale<float>(m_center_offset.x()), - unscale<float>(m_center_offset.y()), 0);
 	        // perform actual slicing
-	        TriangleMeshSlicer mslicer;
 	        const Print *print = this->print();
 	        auto callback = TriangleMeshSlicer::throw_on_cancel_callback_type([print](){print->throw_if_canceled();});
 	        // TriangleMeshSlicer needs the shared vertices.
 	        mesh.require_shared_vertices();
-	        mslicer.init(&mesh, callback);
-	        mslicer.slice(z, mode, float(m_config.slice_closing_radius.value), &layers, callback);
+            MeshSlicingParamsExtended params;
+            params.mode = mode;
+            params.closing_radius = float(m_config.slice_closing_radius.value);
+	        slice_mesh(mesh, z, params, layers, callback);
 	        m_print->throw_if_canceled();
 	    }
 	}
diff --git a/src/libslic3r/SLA/Hollowing.cpp b/src/libslic3r/SLA/Hollowing.cpp
index b38784521..f55a10178 100644
--- a/src/libslic3r/SLA/Hollowing.cpp
+++ b/src/libslic3r/SLA/Hollowing.cpp
@@ -2,6 +2,7 @@
 
 #include <libslic3r/OpenVDBUtils.hpp>
 #include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
 #include <libslic3r/SLA/Hollowing.hpp>
 #include <libslic3r/SLA/IndexedMesh.hpp>
 #include <libslic3r/ClipperUtils.hpp>
@@ -296,10 +297,8 @@ void cut_drainholes(std::vector<ExPolygons> & obj_slices,
     
     mesh.require_shared_vertices();
     
-    TriangleMeshSlicer slicer(&mesh);
-    
     std::vector<ExPolygons> hole_slices;
-    slicer.slice(slicegrid, SlicingMode::Regular, closing_radius, &hole_slices, thr);
+    slice_mesh(mesh, slicegrid, closing_radius, hole_slices, thr);
     
     if (obj_slices.size() != hole_slices.size())
         BOOST_LOG_TRIVIAL(warning)
diff --git a/src/libslic3r/SLA/Pad.cpp b/src/libslic3r/SLA/Pad.cpp
index e11914a1c..658ecf6d8 100644
--- a/src/libslic3r/SLA/Pad.cpp
+++ b/src/libslic3r/SLA/Pad.cpp
@@ -2,6 +2,7 @@
 #include <libslic3r/SLA/SpatIndex.hpp>
 #include <libslic3r/SLA/BoostAdapter.hpp>
 #include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
 
 #include "ConcaveHull.hpp"
 
@@ -476,10 +477,9 @@ void pad_blueprint(const TriangleMesh &      mesh,
                    ThrowOnCancel             thrfn)
 {
     if (mesh.empty()) return;
-    TriangleMeshSlicer slicer(&mesh);
 
     auto out = reserve_vector<ExPolygons>(heights.size());
-    slicer.slice(heights, SlicingMode::Regular, 0.f, &out, thrfn);
+    slice_mesh(mesh, heights, out, thrfn);
 
     size_t count = 0;
     for(auto& o : out) count += o.size();
diff --git a/src/libslic3r/SLA/RasterBase.hpp b/src/libslic3r/SLA/RasterBase.hpp
index bbb83b5a0..1cd688ccb 100644
--- a/src/libslic3r/SLA/RasterBase.hpp
+++ b/src/libslic3r/SLA/RasterBase.hpp
@@ -72,7 +72,8 @@ public:
         size_t width_px = 0;
         size_t height_px = 0;
         
-        Resolution(size_t w = 0, size_t h = 0) : width_px(w), height_px(h) {}
+        Resolution() = default;
+        Resolution(size_t w, size_t h) : width_px(w), height_px(h) {}
         size_t pixels() const { return width_px * height_px; }
     };
     
@@ -81,7 +82,8 @@ public:
         double w_mm = 1.;
         double h_mm = 1.;
         
-        PixelDim(double px_width_mm = 0.0, double px_height_mm = 0.0)
+        PixelDim() = default;
+        PixelDim(double px_width_mm, double px_height_mm)
             : w_mm(px_width_mm), h_mm(px_height_mm)
         {}
     };
diff --git a/src/libslic3r/SLA/SupportTree.cpp b/src/libslic3r/SLA/SupportTree.cpp
index 1bb4cfab7..45e90a704 100644
--- a/src/libslic3r/SLA/SupportTree.cpp
+++ b/src/libslic3r/SLA/SupportTree.cpp
@@ -12,6 +12,7 @@
 #include <libslic3r/MTUtils.hpp>
 #include <libslic3r/ClipperUtils.hpp>
 #include <libslic3r/Model.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
 
 #include <libnest2d/optimizers/nlopt/genetic.hpp>
 #include <libnest2d/optimizers/nlopt/subplex.hpp>
@@ -44,9 +45,7 @@ std::vector<ExPolygons> SupportTree::slice(
 
     if (!sup_mesh.empty()) {
         slices.emplace_back();
-
-        TriangleMeshSlicer sup_slicer(&sup_mesh);
-        sup_slicer.slice(grid, SlicingMode::Regular, cr, &slices.back(), ctl().cancelfn);
+        slice_mesh(sup_mesh, grid, cr, slices.back(), ctl().cancelfn);
     }
 
     if (!pad_mesh.empty()) {
@@ -59,8 +58,7 @@ std::vector<ExPolygons> SupportTree::slice(
         auto padgrid = reserve_vector<float>(size_t(cap > 0 ? cap : 0));
         std::copy(grid.begin(), maxzit, std::back_inserter(padgrid));
 
-        TriangleMeshSlicer pad_slicer(&pad_mesh);
-        pad_slicer.slice(padgrid, SlicingMode::Regular, cr, &slices.back(), ctl().cancelfn);
+        slice_mesh(pad_mesh, padgrid, cr, slices.back(), ctl().cancelfn);
     }
 
     size_t len = grid.size();
diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp
index 108159b89..d2c3c06d2 100644
--- a/src/libslic3r/SLAPrintSteps.cpp
+++ b/src/libslic3r/SLAPrintSteps.cpp
@@ -3,6 +3,7 @@
 #include <libslic3r/Exception.hpp>
 #include <libslic3r/SLAPrintSteps.hpp>
 #include <libslic3r/MeshBoolean.hpp>
+#include <libslic3r/TriangleMeshSlicer.hpp>
 
 // Need the cylinder method for the the drainholes in hollowing step
 #include <libslic3r/SLA/SupportTreeBuilder.hpp>
@@ -198,7 +199,7 @@ static std::vector<bool> create_exclude_mask(
     std::vector<bool> exclude_mask(its.indices.size(), false);
 
     std::vector< std::vector<size_t> > neighbor_index =
-            create_neighbor_index(its);
+            create_vertex_faces_index(its);
 
     auto exclude_neighbors = [&neighbor_index, &exclude_mask](const Vec3i &face)
     {
@@ -470,13 +471,11 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po)
     for(auto it = slindex_it; it != po.m_slice_index.end(); ++it)
         po.m_model_height_levels.emplace_back(it->slice_level());
 
-    TriangleMeshSlicer slicer(&mesh);
-
     po.m_model_slices.clear();
     float closing_r  = float(po.config().slice_closing_radius.value);
     auto  thr        = [this]() { m_print->throw_if_canceled(); };
     auto &slice_grid = po.m_model_height_levels;
-    slicer.slice(slice_grid, SlicingMode::Regular, closing_r, &po.m_model_slices, thr);
+    slice_mesh(mesh, slice_grid, closing_r, po.m_model_slices, thr);
 
     sla::Interior *interior = po.m_hollowing_data ?
                                   po.m_hollowing_data->interior.get() :
@@ -486,9 +485,8 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po)
         TriangleMesh interiormesh = sla::get_mesh(*interior);
         interiormesh.repaired = false;
         interiormesh.repair(true);
-        TriangleMeshSlicer interior_slicer(&interiormesh);
         std::vector<ExPolygons> interior_slices;
-        interior_slicer.slice(slice_grid, SlicingMode::Regular, closing_r, &interior_slices, thr);
+        slice_mesh(interiormesh, slice_grid, closing_r, interior_slices, thr);
 
         sla::ccr::for_each(size_t(0), interior_slices.size(),
                            [&po, &interior_slices] (size_t i) {
diff --git a/src/libslic3r/SVG.cpp b/src/libslic3r/SVG.cpp
index 6bc334eec..6a743b1eb 100644
--- a/src/libslic3r/SVG.cpp
+++ b/src/libslic3r/SVG.cpp
@@ -83,12 +83,6 @@ void SVG::draw(const Lines &lines, std::string stroke, coordf_t stroke_width)
         this->draw(l, stroke, stroke_width);
 }
 
-void SVG::draw(const IntersectionLines &lines, std::string stroke)
-{
-    for (const IntersectionLine &il : lines)
-        this->draw((Line)il, stroke);
-}
-
 void SVG::draw(const ExPolygon &expolygon, std::string fill, const float fill_opacity)
 {
     this->fill = fill;
diff --git a/src/libslic3r/SVG.hpp b/src/libslic3r/SVG.hpp
index da9dca093..314fae9b9 100644
--- a/src/libslic3r/SVG.hpp
+++ b/src/libslic3r/SVG.hpp
@@ -43,8 +43,7 @@ public:
     void draw(const Line &line, std::string stroke = "black", coordf_t stroke_width = 0);
     void draw(const ThickLine &line, const std::string &fill, const std::string &stroke, coordf_t stroke_width = 0);
     void draw(const Lines &lines, std::string stroke = "black", coordf_t stroke_width = 0);
-    void draw(const IntersectionLines &lines, std::string stroke = "black");
-
+    
     void draw(const ExPolygon &expolygon, std::string fill = "grey", const float fill_opacity=1.f);
     void draw_outline(const ExPolygon &polygon, std::string stroke_outer = "black", std::string stroke_holes = "blue", coordf_t stroke_width = 0);
     void draw(const ExPolygons &expolygons, std::string fill = "grey", const float fill_opacity=1.f);
diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp
index 7b5c40466..9ce2d3a7f 100644
--- a/src/libslic3r/TriangleMesh.cpp
+++ b/src/libslic3r/TriangleMesh.cpp
@@ -1,46 +1,28 @@
 #include "Exception.hpp"
 #include "TriangleMesh.hpp"
+#include "TriangleMeshSlicer.hpp"
 #include "ClipperUtils.hpp"
 #include "Geometry.hpp"
-#include "Tesselate.hpp"
+
 #include <libqhullcpp/Qhull.h>
 #include <libqhullcpp/QhullFacetList.h>
 #include <libqhullcpp/QhullVertexSet.h>
+
 #include <cmath>
 #include <deque>
 #include <queue>
-#include <set>
 #include <vector>
-#include <map>
 #include <utility>
 #include <algorithm>
-#include <math.h>
 #include <type_traits>
 
 #include <boost/log/trivial.hpp>
 
-#include <tbb/parallel_for.h>
-
 #include <Eigen/Core>
 #include <Eigen/Dense>
 
-// for SLIC3R_DEBUG_SLICE_PROCESSING
-#include "libslic3r.h"
-
-#if 0
-    #define DEBUG
-    #define _DEBUG
-    #undef NDEBUG
-    #define SLIC3R_DEBUG
-// #define SLIC3R_TRIANGLEMESH_DEBUG
-#endif
-
 #include <assert.h>
 
-#if defined(SLIC3R_DEBUG) || defined(SLIC3R_DEBUG_SLICE_PROCESSING)
-#include "SVG.hpp"
-#endif
-
 namespace Slic3r {
 
 TriangleMesh::TriangleMesh(const Pointf3s &points, const std::vector<Vec3i> &facets) : repaired(false)
@@ -53,7 +35,7 @@ TriangleMesh::TriangleMesh(const Pointf3s &points, const std::vector<Vec3i> &fac
     stl.stats.original_num_facets = stl.stats.number_of_facets;
     stl_allocate(&stl);
 
-	for (uint32_t i = 0; i < stl.stats.number_of_facets; ++ i) {
+    for (uint32_t i = 0; i < stl.stats.number_of_facets; ++ i) {
         stl_facet facet;
         facet.vertex[0] = points[facets[i](0)].cast<float>();
         facet.vertex[1] = points[facets[i](1)].cast<float>();
@@ -104,23 +86,23 @@ TriangleMesh::TriangleMesh(const indexed_triangle_set &M) : repaired(false)
 void TriangleMesh::repair(bool update_shared_vertices)
 {
     if (this->repaired) {
-    	if (update_shared_vertices)
-    		this->require_shared_vertices();
-    	return;
+        if (update_shared_vertices)
+            this->require_shared_vertices();
+        return;
     }
 
     // admesh fails when repairing empty meshes
     if (this->stl.stats.number_of_facets == 0)
-    	return;
+        return;
 
     BOOST_LOG_TRIVIAL(debug) << "TriangleMesh::repair() started";
 
     // checking exact
 #ifdef SLIC3R_TRACE_REPAIR
-	BOOST_LOG_TRIVIAL(trace) << "\tstl_check_faces_exact";
+    BOOST_LOG_TRIVIAL(trace) << "\tstl_check_faces_exact";
 #endif /* SLIC3R_TRACE_REPAIR */
-	assert(stl_validate(&this->stl));
-	stl_check_facets_exact(&stl);
+    assert(stl_validate(&this->stl));
+    stl_check_facets_exact(&stl);
     assert(stl_validate(&this->stl));
     stl.stats.facets_w_1_bad_edge = (stl.stats.connected_facets_2_edge - stl.stats.connected_facets_3_edge);
     stl.stats.facets_w_2_bad_edge = (stl.stats.connected_facets_1_edge - stl.stats.connected_facets_2_edge);
@@ -128,17 +110,17 @@ void TriangleMesh::repair(bool update_shared_vertices)
     
     // checking nearby
     //int last_edges_fixed = 0;
-	float tolerance = (float)stl.stats.shortest_edge;
-	float increment = (float)stl.stats.bounding_diameter / 10000.0f;
+    float tolerance = (float)stl.stats.shortest_edge;
+    float increment = (float)stl.stats.bounding_diameter / 10000.0f;
     int iterations = 2;
     if (stl.stats.connected_facets_3_edge < (int)stl.stats.number_of_facets) {
         for (int i = 0; i < iterations; i++) {
             if (stl.stats.connected_facets_3_edge < (int)stl.stats.number_of_facets) {
                 //printf("Checking nearby. Tolerance= %f Iteration=%d of %d...", tolerance, i + 1, iterations);
 #ifdef SLIC3R_TRACE_REPAIR
-				BOOST_LOG_TRIVIAL(trace) << "\tstl_check_faces_nearby";
+                BOOST_LOG_TRIVIAL(trace) << "\tstl_check_faces_nearby";
 #endif /* SLIC3R_TRACE_REPAIR */
-				stl_check_facets_nearby(&stl, tolerance);
+                stl_check_facets_nearby(&stl, tolerance);
                 //printf("  Fixed %d edges.\n", stl.stats.edges_fixed - last_edges_fixed);
                 //last_edges_fixed = stl.stats.edges_fixed;
                 tolerance += increment;
@@ -155,7 +137,7 @@ void TriangleMesh::repair(bool update_shared_vertices)
         BOOST_LOG_TRIVIAL(trace) << "\tstl_remove_unconnected_facets";
 #endif /* SLIC3R_TRACE_REPAIR */
         stl_remove_unconnected_facets(&stl);
-	    assert(stl_validate(&this->stl));
+        assert(stl_validate(&this->stl));
     }
     
     // fill_holes
@@ -207,7 +189,7 @@ void TriangleMesh::repair(bool update_shared_vertices)
     // and it is risky to generate such a structure once the meshes are shared. Do it now.
     this->its.clear();
     if (update_shared_vertices)
-    	this->require_shared_vertices();
+        this->require_shared_vertices();
 }
 
 float TriangleMesh::volume()
@@ -273,18 +255,18 @@ void TriangleMesh::WriteOBJFile(const char* output_file) const
 void TriangleMesh::scale(float factor)
 {
     stl_scale(&(this->stl), factor);
-	for (stl_vertex& v : this->its.vertices)
-		v *= factor;
+    for (stl_vertex& v : this->its.vertices)
+        v *= factor;
 }
 
 void TriangleMesh::scale(const Vec3d &versor)
 {
     stl_scale_versor(&this->stl, versor.cast<float>());
-	for (stl_vertex& v : this->its.vertices) {
-		v.x() *= versor.x();
-		v.y() *= versor.y();
-		v.z() *= versor.z();
-	}
+    for (stl_vertex& v : this->its.vertices) {
+        v.x() *= versor.x();
+        v.y() *= versor.y();
+        v.z() *= versor.z();
+    }
 }
 
 void TriangleMesh::translate(float x, float y, float z)
@@ -292,9 +274,9 @@ void TriangleMesh::translate(float x, float y, float z)
     if (x == 0.f && y == 0.f && z == 0.f)
         return;
     stl_translate_relative(&(this->stl), x, y, z);
-	stl_vertex shift(x, y, z);
-	for (stl_vertex& v : this->its.vertices)
-		v += shift;
+    stl_vertex shift(x, y, z);
+    for (stl_vertex& v : this->its.vertices)
+        v += shift;
 }
 
 void TriangleMesh::translate(const Vec3f &displacement)
@@ -339,15 +321,15 @@ void TriangleMesh::mirror(const Axis &axis)
     if (axis == X) {
         stl_mirror_yz(&this->stl);
         for (stl_vertex &v : this->its.vertices)
-      		v(0) *= -1.0;
+            v(0) *= -1.0;
     } else if (axis == Y) {
         stl_mirror_xz(&this->stl);
         for (stl_vertex &v : this->its.vertices)
-      		v(1) *= -1.0;
+            v(1) *= -1.0;
     } else if (axis == Z) {
         stl_mirror_xy(&this->stl);
         for (stl_vertex &v : this->its.vertices)
-      		v(2) *= -1.0;
+            v(2) *= -1.0;
     }
 }
 
@@ -355,8 +337,8 @@ void TriangleMesh::transform(const Transform3d& t, bool fix_left_handed)
 {
     stl_transform(&stl, t);
     its_transform(its, t);
-	if (fix_left_handed && t.matrix().block(0, 0, 3, 3).determinant() < 0.) {
-		// Left handed transformation is being applied. It is a good idea to flip the faces and their normals.
+    if (fix_left_handed && t.matrix().block(0, 0, 3, 3).determinant() < 0.) {
+        // Left handed transformation is being applied. It is a good idea to flip the faces and their normals.
         // As for the assert: the repair function would fix the normals, reversing would
         // break them again. The caller should provide a mesh that does not need repair.
         // The repair call is left here so things don't break more than they were.
@@ -365,7 +347,7 @@ void TriangleMesh::transform(const Transform3d& t, bool fix_left_handed)
         stl_reverse_all_facets(&stl);
         this->its.clear();
         this->require_shared_vertices();
-	}
+    }
 }
 
 void TriangleMesh::transform(const Matrix3d& m, bool fix_left_handed)
@@ -522,7 +504,7 @@ ExPolygons TriangleMesh::horizontal_projection() const
     auto                delta = scaled<float>(0.01);
     std::vector<float>  deltas { delta, delta, delta };
     paths.reserve(this->stl.stats.number_of_facets);
-	for (const stl_facet &facet : this->stl.facet_start) {
+    for (const stl_facet &facet : this->stl.facet_start) {
         p.points[0] = Point::new_scale(facet.vertex[0](0), facet.vertex[0](1));
         p.points[1] = Point::new_scale(facet.vertex[1](0), facet.vertex[1](1));
         p.points[2] = Point::new_scale(facet.vertex[2](0), facet.vertex[2](1));
@@ -560,12 +542,12 @@ BoundingBoxf3 TriangleMesh::transformed_bounding_box(const Transform3d &trafo) c
     BoundingBoxf3 bbox;
     if (this->its.vertices.empty()) {
         // Using the STL faces.
-		for (const stl_facet &facet : this->stl.facet_start)
+        for (const stl_facet &facet : this->stl.facet_start)
             for (size_t j = 0; j < 3; ++ j)
                 bbox.merge(trafo * facet.vertex[j].cast<double>());
     } else {
         // Using the shared vertices should be a bit quicker than using the STL faces.
-		for (const stl_vertex &v : this->its.vertices)
+        for (const stl_vertex &v : this->its.vertices)
             bbox.merge(trafo * v.cast<double>());
     }
     return bbox;
@@ -576,29 +558,29 @@ TriangleMesh TriangleMesh::convex_hull_3d() const
     // The qhull call:
     orgQhull::Qhull qhull;
     qhull.disableOutputStream(); // we want qhull to be quiet
-	std::vector<realT> src_vertices;
-	try
+    std::vector<realT> src_vertices;
+    try
     {
-    	if (this->has_shared_vertices()) {
+        if (this->has_shared_vertices()) {
 #if REALfloat
-	    	qhull.runQhull("", 3, (int)this->its.vertices.size(), (const realT*)(this->its.vertices.front().data()), "Qt");
+            qhull.runQhull("", 3, (int)this->its.vertices.size(), (const realT*)(this->its.vertices.front().data()), "Qt");
 #else
-	    	src_vertices.reserve(this->its.vertices.size() * 3);
-	    	// We will now fill the vector with input points for computation:
-			for (const stl_vertex &v : this->its.vertices)
-				for (int i = 0; i < 3; ++ i)
-		        	src_vertices.emplace_back(v(i));
-	        qhull.runQhull("", 3, (int)src_vertices.size() / 3, src_vertices.data(), "Qt");
+            src_vertices.reserve(this->its.vertices.size() * 3);
+            // We will now fill the vector with input points for computation:
+            for (const stl_vertex &v : this->its.vertices)
+                for (int i = 0; i < 3; ++ i)
+                    src_vertices.emplace_back(v(i));
+            qhull.runQhull("", 3, (int)src_vertices.size() / 3, src_vertices.data(), "Qt");
 #endif
-	    } else {
-	    	src_vertices.reserve(this->stl.facet_start.size() * 9);
-	    	// We will now fill the vector with input points for computation:
-			for (const stl_facet &f : this->stl.facet_start)
-				for (int i = 0; i < 3; ++ i)
-					for (int j = 0; j < 3; ++ j)
-		        		src_vertices.emplace_back(f.vertex[i](j));
-	        qhull.runQhull("", 3, (int)src_vertices.size() / 3, src_vertices.data(), "Qt");
-	    }
+        } else {
+            src_vertices.reserve(this->stl.facet_start.size() * 9);
+            // We will now fill the vector with input points for computation:
+            for (const stl_facet &f : this->stl.facet_start)
+                for (int i = 0; i < 3; ++ i)
+                    for (int j = 0; j < 3; ++ j)
+                        src_vertices.emplace_back(f.vertex[i](j));
+            qhull.runQhull("", 3, (int)src_vertices.size() / 3, src_vertices.data(), "Qt");
+        }
     }
     catch (...)
     {
@@ -633,9 +615,8 @@ std::vector<ExPolygons> TriangleMesh::slice(const std::vector<double> &z)
 {
     // convert doubles to floats
     std::vector<float> z_f(z.begin(), z.end());
-    TriangleMeshSlicer mslicer(this);
     std::vector<ExPolygons> layers;
-    mslicer.slice(z_f, SlicingMode::Regular, 0.0004f, &layers, [](){});
+    slice_mesh(*this, z_f, 0.0004f, layers);
     return layers;
 }
 
@@ -655,49 +636,61 @@ void TriangleMesh::require_shared_vertices()
 
 size_t TriangleMesh::memsize() const
 {
-	size_t memsize = 8 + this->stl.memsize() + this->its.memsize();
-	return memsize;
+    size_t memsize = 8 + this->stl.memsize() + this->its.memsize();
+    return memsize;
 }
 
 // Release optional data from the mesh if the object is on the Undo / Redo stack only. Returns the amount of memory released.
 size_t TriangleMesh::release_optional()
 {
-	size_t memsize_released = sizeof(stl_neighbors) * this->stl.neighbors_start.size() + this->its.memsize();
-	// The indexed triangle set may be recalculated using the stl_generate_shared_vertices() function.
-	this->its.clear();
-	// The neighbors structure may be recalculated using the stl_check_facets_exact() function.
-	this->stl.neighbors_start.clear();
-	return memsize_released;
+    size_t memsize_released = sizeof(stl_neighbors) * this->stl.neighbors_start.size() + this->its.memsize();
+    // The indexed triangle set may be recalculated using the stl_generate_shared_vertices() function.
+    this->its.clear();
+    // The neighbors structure may be recalculated using the stl_check_facets_exact() function.
+    this->stl.neighbors_start.clear();
+    return memsize_released;
 }
 
 // Restore optional data possibly released by release_optional().
 void TriangleMesh::restore_optional()
 {
-	if (! this->stl.facet_start.empty()) {
-		// Save the old stats before calling stl_check_faces_exact, as it may modify the statistics.
-		stl_stats stats = this->stl.stats;
-		if (this->stl.neighbors_start.empty()) {
-			stl_reallocate(&this->stl);
-			stl_check_facets_exact(&this->stl);
-		}
-		if (this->its.vertices.empty())
-			stl_generate_shared_vertices(&this->stl, this->its);
-		// Restore the old statistics.
-		this->stl.stats = stats;
-	}
+    if (! this->stl.facet_start.empty()) {
+        // Save the old stats before calling stl_check_faces_exact, as it may modify the statistics.
+        stl_stats stats = this->stl.stats;
+        if (this->stl.neighbors_start.empty()) {
+            stl_reallocate(&this->stl);
+            stl_check_facets_exact(&this->stl);
+        }
+        if (this->its.vertices.empty())
+            stl_generate_shared_vertices(&this->stl, this->its);
+        // Restore the old statistics.
+        this->stl.stats = stats;
+    }
 }
 
-void TriangleMeshSlicer::init(const TriangleMesh *_mesh, throw_on_cancel_callback_type throw_on_cancel)
+std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangle_set &its)
 {
-    mesh = _mesh;
-    if (! mesh->has_shared_vertices())
-        throw Slic3r::InvalidArgument("TriangleMeshSlicer was passed a mesh without shared vertices.");
+    std::vector<std::vector<size_t>> index;
 
-    throw_on_cancel();
-    facets_edges.assign(_mesh->stl.stats.number_of_facets * 3, -1);
-	v_scaled_shared.assign(_mesh->its.vertices.size(), stl_vertex());
-	for (size_t i = 0; i < v_scaled_shared.size(); ++ i)
-        this->v_scaled_shared[i] = _mesh->its.vertices[i] / float(SCALING_FACTOR);
+    if (! its.vertices.empty()) {
+        size_t res = its.indices.size() / its.vertices.size();
+        index.assign(its.vertices.size(), reserve_vector<size_t>(res));
+        for (size_t fi = 0; fi < its.indices.size(); ++fi) {
+            auto &face = its.indices[fi];
+            index[face(0)].emplace_back(fi);
+            index[face(1)].emplace_back(fi);
+            index[face(2)].emplace_back(fi);
+        }
+    }
+
+    return index;
+}
+
+// Map from a facet edge to a neighbor face index or -1 if no neighbor exists.
+template<typename ThrowOnCancelCallback>
+static inline std::vector<int> create_face_neighbors_index_impl(const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
+{
+    std::vector<int> out(its.indices.size() * 3, -1);
 
     // Create a mapping from triangle edge into face.
     struct EdgeToFace {
@@ -713,12 +706,12 @@ void TriangleMeshSlicer::init(const TriangleMesh *_mesh, throw_on_cancel_callbac
         bool operator<(const EdgeToFace &other) const { return vertex_low < other.vertex_low || (vertex_low == other.vertex_low && vertex_high < other.vertex_high); }
     };
     std::vector<EdgeToFace> edges_map;
-    edges_map.assign(this->mesh->stl.stats.number_of_facets * 3, EdgeToFace());
-    for (uint32_t facet_idx = 0; facet_idx < this->mesh->stl.stats.number_of_facets; ++ facet_idx)
+    edges_map.assign(its.indices.size() * 3, EdgeToFace());
+    for (uint32_t facet_idx = 0; facet_idx < its.indices.size(); ++ facet_idx)
         for (int i = 0; i < 3; ++ i) {
-            EdgeToFace &e2f = edges_map[facet_idx*3+i];
-            e2f.vertex_low  = this->mesh->its.indices[facet_idx][i];
-            e2f.vertex_high = this->mesh->its.indices[facet_idx][(i + 1) % 3];
+            EdgeToFace &e2f = edges_map[facet_idx * 3 + i];
+            e2f.vertex_low  = its.indices[facet_idx][i];
+            e2f.vertex_high = its.indices[facet_idx][(i + 1) % 3];
             e2f.face        = facet_idx;
             // 1 based indexing, to be always strictly positive.
             e2f.face_edge   = i + 1;
@@ -761,10 +754,10 @@ void TriangleMeshSlicer::init(const TriangleMesh *_mesh, throw_on_cancel_callbac
                 }
         }
         // Assign an edge index to the 1st face.
-        this->facets_edges[edge_i.face * 3 + std::abs(edge_i.face_edge) - 1] = num_edges;
+        out[edge_i.face * 3 + std::abs(edge_i.face_edge) - 1] = num_edges;
         if (found) {
             EdgeToFace &edge_j = edges_map[j];
-            this->facets_edges[edge_j.face * 3 + std::abs(edge_j.face_edge) - 1] = num_edges;
+            out[edge_j.face * 3 + std::abs(edge_j.face_edge) - 1] = num_edges;
             // Mark the edge as connected.
             edge_j.face = -1;
         }
@@ -772,1179 +765,75 @@ void TriangleMeshSlicer::init(const TriangleMesh *_mesh, throw_on_cancel_callbac
         if ((i & 0x0ffff) == 0)
             throw_on_cancel();
     }
-}
 
-
-
-void TriangleMeshSlicer::set_up_direction(const Vec3f& up)
-{
-    m_quaternion.setFromTwoVectors(up, Vec3f::UnitZ());
-    m_use_quaternion = true;
-}
-
-void TriangleMeshSlicer::slice(
-    const std::vector<float> &z, 
-    SlicingMode mode, size_t alternate_mode_first_n_layers, SlicingMode alternate_mode,
-    std::vector<Polygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const
-{
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice";
-
-    /*
-       This method gets called with a list of unscaled Z coordinates and outputs
-       a vector pointer having the same number of items as the original list.
-       Each item is a vector of polygons created by slicing our mesh at the 
-       given heights.
-       
-       This method should basically combine the behavior of the existing
-       Perl methods defined in lib/Slic3r/TriangleMesh.pm:
-       
-       - analyze(): this creates the 'facets_edges' and the 'edges_facets'
-            tables (we don't need the 'edges' table)
-       
-       - slice_facet(): this has to be done for each facet. It generates 
-            intersection lines with each plane identified by the Z list.
-            The get_layer_range() binary search used to identify the Z range
-            of the facet is already ported to C++ (see Object.xsp)
-       
-       - make_loops(): this has to be done for each layer. It creates polygons
-            from the lines generated by the previous step.
-        
-        At the end, we free the tables generated by analyze() as we don't 
-        need them anymore.
-        
-        NOTE: this method accepts a vector of floats because the mesh coordinate
-        type is float.
-    */
-    
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::_slice_do";
-    std::vector<IntersectionLines> lines(z.size());
-    {
-        boost::mutex lines_mutex;
-        tbb::parallel_for(
-            tbb::blocked_range<int>(0,this->mesh->stl.stats.number_of_facets),
-            [&lines, &lines_mutex, &z, throw_on_cancel, this](const tbb::blocked_range<int>& range) {
-                for (int facet_idx = range.begin(); facet_idx < range.end(); ++ facet_idx) {
-                    if ((facet_idx & 0x0ffff) == 0)
-                        throw_on_cancel();
-                    this->_slice_do(facet_idx, &lines, &lines_mutex, z);
-                }
-            }
-        );
-    }
-    throw_on_cancel();
-
-    // v_scaled_shared could be freed here
-    
-    // build loops
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::_make_loops_do";
-    layers->resize(z.size());
-    tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, z.size()),
-        [&lines, &layers, mode, alternate_mode_first_n_layers, alternate_mode, throw_on_cancel, this](const tbb::blocked_range<size_t>& range) {
-            for (size_t line_idx = range.begin(); line_idx < range.end(); ++ line_idx) {
-                if ((line_idx & 0x0ffff) == 0)
-                    throw_on_cancel();
-
-                Polygons   &polygons  = (*layers)[line_idx];
-                this->make_loops(lines[line_idx], &polygons);
-
-                auto this_mode = line_idx < alternate_mode_first_n_layers ? alternate_mode : mode;
-                if (! polygons.empty()) {
-                    if (this_mode == SlicingMode::Positive) {
-                        // Reorient all loops to be CCW.
-                        for (Polygon& p : polygons)
-                            p.make_counter_clockwise();
-                    } else if (this_mode == SlicingMode::PositiveLargestContour) {
-                        // Keep just the largest polygon, make it CCW.
-                        double   max_area = 0.;
-                        Polygon* max_area_polygon = nullptr;
-                        for (Polygon& p : polygons) {
-                            double a = p.area();
-                            if (std::abs(a) > std::abs(max_area)) {
-                                max_area = a;
-                                max_area_polygon = &p;
-                            }
-                        }
-                        assert(max_area_polygon != nullptr);
-                        if (max_area < 0.)
-                            max_area_polygon->reverse();
-                        Polygon p(std::move(*max_area_polygon));
-                        polygons.clear();
-                        polygons.emplace_back(std::move(p));
-                    }
-                }
-            }
-        }
-    );
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice finished";
-
-#ifdef SLIC3R_DEBUG
-    {
-        static int iRun = 0;
-        for (size_t i = 0; i < z.size(); ++ i) {
-            Polygons  &polygons   = (*layers)[i];
-            ExPolygons expolygons = union_ex(polygons, true);
-            SVG::export_expolygons(debug_out_path("slice_%d_%d.svg", iRun, i).c_str(), expolygons);
-            {
-                BoundingBox bbox;
-                for (const IntersectionLine &l : lines[i]) {
-                    bbox.merge(l.a);
-                    bbox.merge(l.b);
-                }
-                SVG svg(debug_out_path("slice_loops_%d_%d.svg", iRun, i).c_str(), bbox);
-                svg.draw(expolygons);
-                for (const IntersectionLine &l : lines[i])
-                    svg.draw(l, "red", 0);
-                svg.draw_outline(expolygons, "black", "blue", 0);
-                svg.Close();
-            }
-#if 0
-//FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
-            for (Polygon &poly : polygons) {
-                for (size_t i = 1; i < poly.points.size(); ++ i)
-                    assert(poly.points[i-1] != poly.points[i]);
-                assert(poly.points.front() != poly.points.back());
-            }
-#endif
-        }
-        ++ iRun;
-    }
-#endif
-}
-
-void TriangleMeshSlicer::_slice_do(size_t facet_idx, std::vector<IntersectionLines>* lines, boost::mutex* lines_mutex, 
-    const std::vector<float> &z) const
-{
-    const stl_facet &facet = m_use_quaternion ? (this->mesh->stl.facet_start.data() + facet_idx)->rotated(m_quaternion) : *(this->mesh->stl.facet_start.data() + facet_idx);
-    
-    // find facet extents
-    const float min_z = fminf(facet.vertex[0](2), fminf(facet.vertex[1](2), facet.vertex[2](2)));
-    const float max_z = fmaxf(facet.vertex[0](2), fmaxf(facet.vertex[1](2), facet.vertex[2](2)));
-    
-    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
-    printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx,
-        facet.vertex[0](0), facet.vertex[0](1), facet.vertex[0](2),
-        facet.vertex[1](0), facet.vertex[1](1), facet.vertex[1](2),
-        facet.vertex[2](0), facet.vertex[2](1), facet.vertex[2](2));
-    printf("z: min = %.2f, max = %.2f\n", min_z, max_z);
-    #endif /* SLIC3R_TRIANGLEMESH_DEBUG */
-    
-    // find layer extents
-    std::vector<float>::const_iterator min_layer, max_layer;
-    min_layer = std::lower_bound(z.begin(), z.end(), min_z); // first layer whose slice_z is >= min_z
-    max_layer = std::upper_bound(min_layer, z.end(), max_z); // first layer whose slice_z is > max_z
-    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
-    printf("layers: min = %d, max = %d\n", (int)(min_layer - z.begin()), (int)(max_layer - z.begin()));
-    #endif /* SLIC3R_TRIANGLEMESH_DEBUG */
-    
-    for (std::vector<float>::const_iterator it = min_layer; it != max_layer; ++ it) {
-        std::vector<float>::size_type layer_idx = it - z.begin();
-        IntersectionLine il;
-        if (this->slice_facet(*it / SCALING_FACTOR, facet, facet_idx, min_z, max_z, &il) == TriangleMeshSlicer::Slicing) {
-            boost::lock_guard<boost::mutex> l(*lines_mutex);
-            if (il.edge_type == feHorizontal) {
-                // Ignore horizontal triangles. Any valid horizontal triangle must have a vertical triangle connected, otherwise the part has zero volume.
-            } else
-                (*lines)[layer_idx].emplace_back(il);
-        }
-    }
-}
-
-void TriangleMeshSlicer::slice(
-    const std::vector<float> &z, SlicingMode mode, size_t alternate_mode_first_n_layers, SlicingMode alternate_mode, const float closing_radius, 
-    std::vector<ExPolygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const
-{
-    std::vector<Polygons> layers_p;
-    this->slice(z, 
-        (mode == SlicingMode::PositiveLargestContour) ? SlicingMode::Positive : mode, 
-        alternate_mode_first_n_layers,
-        (alternate_mode == SlicingMode::PositiveLargestContour) ? SlicingMode::Positive : alternate_mode,
-        &layers_p, throw_on_cancel);
-    
-	BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - start";
-	layers->resize(z.size());
-	tbb::parallel_for(
-		tbb::blocked_range<size_t>(0, z.size()),
-		[&layers_p, mode, alternate_mode_first_n_layers, alternate_mode, closing_radius, layers, throw_on_cancel, this]
-        (const tbb::blocked_range<size_t>& range) {
-    		for (size_t layer_id = range.begin(); layer_id < range.end(); ++ layer_id) {
-#ifdef SLIC3R_TRIANGLEMESH_DEBUG
-                printf("Layer %zu (slice_z = %.2f):\n", layer_id, z[layer_id]);
-#endif
-                throw_on_cancel();
-                ExPolygons &expolygons = (*layers)[layer_id];
-    			this->make_expolygons(layers_p[layer_id], closing_radius, &expolygons);
-                const auto this_mode = layer_id < alternate_mode_first_n_layers ? alternate_mode : mode;
-    			if (this_mode == SlicingMode::PositiveLargestContour)
-					keep_largest_contour_only(expolygons);
-    		}
-    	});
-	BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - end";
-}
-
-// Return true, if the facet has been sliced and line_out has been filled.
-TriangleMeshSlicer::FacetSliceType TriangleMeshSlicer::slice_facet(
-    float slice_z, const stl_facet &facet, const int facet_idx,
-    const float min_z, const float max_z, 
-    IntersectionLine *line_out) const
-{
-    IntersectionPoint points[3];
-    size_t            num_points = 0;
-    size_t            point_on_layer = size_t(-1);
-
-    // Reorder vertices so that the first one is the one with lowest Z.
-    // This is needed to get all intersection lines in a consistent order
-    // (external on the right of the line)
-    const stl_triangle_vertex_indices &vertices = this->mesh->its.indices[facet_idx];
-    int i = (facet.vertex[1].z() == min_z) ? 1 : ((facet.vertex[2].z() == min_z) ? 2 : 0);
-
-    // These are used only if the cut plane is tilted:
-    stl_vertex rotated_a;
-    stl_vertex rotated_b;
-
-    for (int j = i; j - i < 3; ++j) {  // loop through facet edges
-        int        edge_id  = this->facets_edges[facet_idx * 3 + (j % 3)];
-        int        a_id     = vertices[j % 3];
-        int        b_id     = vertices[(j+1) % 3];
-
-        const stl_vertex *a;
-        const stl_vertex *b;
-        if (m_use_quaternion) {
-            rotated_a = m_quaternion * this->v_scaled_shared[a_id];
-            rotated_b = m_quaternion * this->v_scaled_shared[b_id];
-            a = &rotated_a;
-            b = &rotated_b;
-        }
-        else {
-            a = &this->v_scaled_shared[a_id];
-            b = &this->v_scaled_shared[b_id];
-        }
-        
-        // Is edge or face aligned with the cutting plane?
-        if (a->z() == slice_z && b->z() == slice_z) {
-            // Edge is horizontal and belongs to the current layer.
-            // The following rotation of the three vertices may not be efficient, but this branch happens rarely.
-            const stl_vertex &v0 = m_use_quaternion ? stl_vertex(m_quaternion * this->v_scaled_shared[vertices[0]]) : this->v_scaled_shared[vertices[0]];
-            const stl_vertex &v1 = m_use_quaternion ? stl_vertex(m_quaternion * this->v_scaled_shared[vertices[1]]) : this->v_scaled_shared[vertices[1]];
-            const stl_vertex &v2 = m_use_quaternion ? stl_vertex(m_quaternion * this->v_scaled_shared[vertices[2]]) : this->v_scaled_shared[vertices[2]];
-            const stl_normal &normal = facet.normal;
-            // We may ignore this edge for slicing purposes, but we may still use it for object cutting.
-            FacetSliceType    result = Slicing;
-            if (min_z == max_z) {
-                // All three vertices are aligned with slice_z.
-                line_out->edge_type = feHorizontal;
-                result = Cutting;
-                if (normal.z() < 0) {
-                    // If normal points downwards this is a bottom horizontal facet so we reverse its point order.
-                    std::swap(a, b);
-                    std::swap(a_id, b_id);
-                }
-            } else {
-                // Two vertices are aligned with the cutting plane, the third vertex is below or above the cutting plane.
-                // Is the third vertex below the cutting plane?
-                bool third_below = v0.z() < slice_z || v1.z() < slice_z || v2.z() < slice_z;
-                // Two vertices on the cutting plane, the third vertex is below the plane. Consider the edge to be part of the slice
-                // only if it is the upper edge.
-                // (the bottom most edge resp. vertex of a triangle is not owned by the triangle, but the top most edge resp. vertex is part of the triangle
-                // in respect to the cutting plane).
-                result = third_below ? Slicing : Cutting;
-                if (third_below) {
-                    line_out->edge_type = feTop;
-                    std::swap(a, b);
-                    std::swap(a_id, b_id);
-                } else
-                    line_out->edge_type = feBottom;
-            }
-            line_out->a.x()  = a->x();
-            line_out->a.y()  = a->y();
-            line_out->b.x()  = b->x();
-            line_out->b.y()  = b->y();
-            line_out->a_id   = a_id;
-            line_out->b_id   = b_id;
-            assert(line_out->a != line_out->b);
-            return result;
-        }
-
-        if (a->z() == slice_z) {
-            // Only point a alings with the cutting plane.
-            if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != a_id) {
-                point_on_layer = num_points;
-                IntersectionPoint &point = points[num_points ++];
-                point.x()      = a->x();
-                point.y()      = a->y();
-                point.point_id = a_id;
-            }
-        } else if (b->z() == slice_z) {
-            // Only point b alings with the cutting plane.
-            if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != b_id) {
-                point_on_layer = num_points;
-                IntersectionPoint &point = points[num_points ++];
-                point.x()      = b->x();
-                point.y()      = b->y();
-                point.point_id = b_id;
-            }
-        } else if ((a->z() < slice_z && b->z() > slice_z) || (b->z() < slice_z && a->z() > slice_z)) {
-            // A general case. The face edge intersects the cutting plane. Calculate the intersection point.
-            assert(a_id != b_id);
-            // Sort the edge to give a consistent answer.
-            if (a_id > b_id) {
-                std::swap(a_id, b_id);
-                std::swap(a, b);
-            }
-            IntersectionPoint &point = points[num_points];
-			double t = (double(slice_z) - double(b->z())) / (double(a->z()) - double(b->z()));
-            if (t <= 0.) {
-                if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != a_id) {
-                    point.x() = a->x();
-                    point.y() = a->y();
-                    point_on_layer = num_points ++;
-                    point.point_id = a_id;
-                }
-            } else if (t >= 1.) {
-                if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != b_id) {
-                    point.x() = b->x();
-                    point.y() = b->y();
-                    point_on_layer = num_points ++;
-                    point.point_id = b_id;
-                }
-            } else {
-                point.x() = coord_t(floor(double(b->x()) + (double(a->x()) - double(b->x())) * t + 0.5));
-                point.y() = coord_t(floor(double(b->y()) + (double(a->y()) - double(b->y())) * t + 0.5));
-                point.edge_id = edge_id;
-                ++ num_points;
-            }
-        }
-    }
-
-    // Facets must intersect each plane 0 or 2 times, or it may touch the plane at a single vertex only.
-    assert(num_points < 3);
-    if (num_points == 2) {
-        line_out->edge_type  = feGeneral;
-        line_out->a          = (Point)points[1];
-        line_out->b          = (Point)points[0];
-        line_out->a_id       = points[1].point_id;
-        line_out->b_id       = points[0].point_id;
-        line_out->edge_a_id  = points[1].edge_id;
-        line_out->edge_b_id  = points[0].edge_id;
-        // Not a zero lenght edge.
-        //FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
-        //assert(line_out->a != line_out->b);
-        // The plane cuts at least one edge in a general position.
-        assert(line_out->a_id == -1 || line_out->b_id == -1);
-        assert(line_out->edge_a_id != -1 || line_out->edge_b_id != -1);
-        // General slicing position, use the segment for both slicing and object cutting.
-#if 0
-        if (line_out->a_id != -1 && line_out->b_id != -1) {
-            // Solving a degenerate case, where both the intersections snapped to an edge.
-            // Correctly classify the face as below or above based on the position of the 3rd point.
-            int i = vertices[0];
-            if (i == line_out->a_id || i == line_out->b_id)
-                i = vertices[1];
-            if (i == line_out->a_id || i == line_out->b_id)
-                i = vertices[2];
-            assert(i != line_out->a_id && i != line_out->b_id);
-            line_out->edge_type = ((m_use_quaternion ?
-                                    (m_quaternion * this->v_scaled_shared[i]).z()
-                                    : this->v_scaled_shared[i].z()) < slice_z) ? feTop : feBottom;
-        }
-#endif
-        return Slicing;
-    }
-    return NoSlice;
-}
-
-#if 0
-//FIXME Should this go away? For valid meshes the function slice_facet() returns Slicing
-// and sets edges of vertical triangles to produce only a single edge per pair of neighbor faces.
-// So the following code makes only sense now to handle degenerate meshes with more than two faces
-// sharing a single edge.
-static inline void remove_tangent_edges(std::vector<IntersectionLine> &lines)
-{
-    std::vector<IntersectionLine*> by_vertex_pair;
-    by_vertex_pair.reserve(lines.size());
-    for (IntersectionLine& line : lines)
-        if (line.edge_type != feGeneral && line.a_id != -1)
-            // This is a face edge. Check whether there is its neighbor stored in lines.
-            by_vertex_pair.emplace_back(&line);
-    auto edges_lower_sorted = [](const IntersectionLine *l1, const IntersectionLine *l2) {
-        // Sort vertices of l1, l2 lexicographically
-        int l1a = l1->a_id;
-        int l1b = l1->b_id;
-        int l2a = l2->a_id;
-        int l2b = l2->b_id;
-        if (l1a > l1b)
-            std::swap(l1a, l1b);
-        if (l2a > l2b)
-            std::swap(l2a, l2b);
-        // Lexicographical "lower" operator on lexicographically sorted vertices should bring equal edges together when sored.
-        return l1a < l2a || (l1a == l2a && l1b < l2b);
-    };
-    std::sort(by_vertex_pair.begin(), by_vertex_pair.end(), edges_lower_sorted);
-    for (auto line = by_vertex_pair.begin(); line != by_vertex_pair.end(); ++ line) {
-        IntersectionLine &l1 = **line;
-        if (! l1.skip()) {
-            // Iterate as long as line and line2 edges share the same end points.
-            for (auto line2 = line + 1; line2 != by_vertex_pair.end() && ! edges_lower_sorted(*line, *line2); ++ line2) {
-                // Lines must share the end points.
-                assert(! edges_lower_sorted(*line, *line2));
-                assert(! edges_lower_sorted(*line2, *line));
-                IntersectionLine &l2 = **line2;
-                if (l2.skip())
-                    continue;
-                if (l1.a_id == l2.a_id) {
-                    assert(l1.b_id == l2.b_id);
-                    l2.set_skip();
-                    // If they are both oriented upwards or downwards (like a 'V'),
-                    // then we can remove both edges from this layer since it won't 
-                    // affect the sliced shape.
-                    // If one of them is oriented upwards and the other is oriented
-                    // downwards, let's only keep one of them (it doesn't matter which
-                    // one since all 'top' lines were reversed at slicing).
-                    if (l1.edge_type == l2.edge_type) {
-                        l1.set_skip();
-                        break;
-                    }
-                } else {
-                    assert(l1.a_id == l2.b_id && l1.b_id == l2.a_id);
-                    // If this edge joins two horizontal facets, remove both of them.
-                    if (l1.edge_type == feHorizontal && l2.edge_type == feHorizontal) {
-                        l1.set_skip();
-                        l2.set_skip();
-                        break;
-                    }
-                }
-            }
-        }
-    }
-}
-#endif
-
-struct OpenPolyline {
-    OpenPolyline() {};
-    OpenPolyline(const IntersectionReference &start, const IntersectionReference &end, Points &&points) : 
-        start(start), end(end), points(std::move(points)), consumed(false) { this->length = Slic3r::length(this->points); }
-    void reverse() {
-        std::swap(start, end);
-        std::reverse(points.begin(), points.end());
-    }
-    IntersectionReference   start;
-    IntersectionReference   end;
-    Points                  points;
-    double                  length;
-    bool                    consumed;
-};
-
-// called by TriangleMeshSlicer::make_loops() to connect sliced triangles into closed loops and open polylines by the triangle connectivity.
-// Only connects segments crossing triangles of the same orientation.
-static void chain_lines_by_triangle_connectivity(std::vector<IntersectionLine> &lines, Polygons &loops, std::vector<OpenPolyline> &open_polylines)
-{
-    // Build a map of lines by edge_a_id and a_id.
-    std::vector<IntersectionLine*> by_edge_a_id;
-    std::vector<IntersectionLine*> by_a_id;
-    by_edge_a_id.reserve(lines.size());
-    by_a_id.reserve(lines.size());
-    for (IntersectionLine &line : lines) {
-        if (! line.skip()) {
-            if (line.edge_a_id != -1)
-                by_edge_a_id.emplace_back(&line);
-            if (line.a_id != -1)
-                by_a_id.emplace_back(&line);
-        }
-    }
-    auto by_edge_lower = [](const IntersectionLine* il1, const IntersectionLine *il2) { return il1->edge_a_id < il2->edge_a_id; };
-    auto by_vertex_lower = [](const IntersectionLine* il1, const IntersectionLine *il2) { return il1->a_id < il2->a_id; };
-    std::sort(by_edge_a_id.begin(), by_edge_a_id.end(), by_edge_lower);
-    std::sort(by_a_id.begin(), by_a_id.end(), by_vertex_lower);
-    // Chain the segments with a greedy algorithm, collect the loops and unclosed polylines.
-    IntersectionLines::iterator it_line_seed = lines.begin();
-    for (;;) {
-        // take first spare line and start a new loop
-        IntersectionLine *first_line = nullptr;
-        for (; it_line_seed != lines.end(); ++ it_line_seed)
-            if (it_line_seed->is_seed_candidate()) {
-            //if (! it_line_seed->skip()) {
-                first_line = &(*it_line_seed ++);
-                break;
-            }
-        if (first_line == nullptr)
-            break;
-        first_line->set_skip();
-        Points loop_pts;
-        loop_pts.emplace_back(first_line->a);
-        IntersectionLine *last_line = first_line;
-        
-        /*
-        printf("first_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", 
-            first_line->edge_a_id, first_line->edge_b_id, first_line->a_id, first_line->b_id,
-            first_line->a.x, first_line->a.y, first_line->b.x, first_line->b.y);
-        */
-        
-        IntersectionLine key;
-        for (;;) {
-            // find a line starting where last one finishes
-            IntersectionLine* next_line = nullptr;
-            if (last_line->edge_b_id != -1) {
-                key.edge_a_id = last_line->edge_b_id;
-                auto it_begin = std::lower_bound(by_edge_a_id.begin(), by_edge_a_id.end(), &key, by_edge_lower);
-                if (it_begin != by_edge_a_id.end()) {
-                    auto it_end = std::upper_bound(it_begin, by_edge_a_id.end(), &key, by_edge_lower);
-                    for (auto it_line = it_begin; it_line != it_end; ++ it_line)
-                        if (! (*it_line)->skip()) {
-                            next_line = *it_line;
-                            break;
-                        }
-                }
-            }
-            if (next_line == nullptr && last_line->b_id != -1) {
-                key.a_id = last_line->b_id;
-                auto it_begin = std::lower_bound(by_a_id.begin(), by_a_id.end(), &key, by_vertex_lower);
-                if (it_begin != by_a_id.end()) {
-                    auto it_end = std::upper_bound(it_begin, by_a_id.end(), &key, by_vertex_lower);
-                    for (auto it_line = it_begin; it_line != it_end; ++ it_line)
-                        if (! (*it_line)->skip()) {
-                            next_line = *it_line;
-                            break;
-                        }
-                }
-            }
-            if (next_line == nullptr) {
-                // Check whether we closed this loop.
-                if ((first_line->edge_a_id != -1 && first_line->edge_a_id == last_line->edge_b_id) || 
-                    (first_line->a_id      != -1 && first_line->a_id      == last_line->b_id)) {
-                    // The current loop is complete. Add it to the output.
-                    loops.emplace_back(std::move(loop_pts));
-                    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
-                    printf("  Discovered %s polygon of %d points\n", (p.is_counter_clockwise() ? "ccw" : "cw"), (int)p.points.size());
-                    #endif
-                } else {
-                    // This is an open polyline. Add it to the list of open polylines. These open polylines will processed later.
-                    loop_pts.emplace_back(last_line->b);
-                    open_polylines.emplace_back(OpenPolyline(
-                        IntersectionReference(first_line->a_id, first_line->edge_a_id), 
-                        IntersectionReference(last_line->b_id, last_line->edge_b_id), std::move(loop_pts)));
-                }
-                break;
-            }
-            /*
-            printf("next_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", 
-                next_line->edge_a_id, next_line->edge_b_id, next_line->a_id, next_line->b_id,
-                next_line->a.x, next_line->a.y, next_line->b.x, next_line->b.y);
-            */
-            loop_pts.emplace_back(next_line->a);
-            last_line = next_line;
-            next_line->set_skip();
-        }
-    }
-}
-
-std::vector<OpenPolyline*> open_polylines_sorted(std::vector<OpenPolyline> &open_polylines, bool update_lengths)
-{
-    std::vector<OpenPolyline*> out;
-    out.reserve(open_polylines.size());
-    for (OpenPolyline &opl : open_polylines)
-        if (! opl.consumed) {
-            if (update_lengths)
-                opl.length = Slic3r::length(opl.points);
-            out.emplace_back(&opl);
-        }
-    std::sort(out.begin(), out.end(), [](const OpenPolyline *lhs, const OpenPolyline *rhs){ return lhs->length > rhs->length; });
     return out;
 }
 
-// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices.
-// Depending on "try_connect_reversed", it may or may not connect segments crossing triangles of opposite orientation.
-static void chain_open_polylines_exact(std::vector<OpenPolyline> &open_polylines, Polygons &loops, bool try_connect_reversed)
+std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its)
 {
-    // Store the end points of open_polylines into vectors sorted
-    struct OpenPolylineEnd {
-        OpenPolylineEnd(OpenPolyline *polyline, bool start) : polyline(polyline), start(start) {}
-        OpenPolyline    *polyline;
-        // Is it the start or end point?
-        bool             start;
-        const IntersectionReference& ipref() const { return start ? polyline->start : polyline->end; }
-        // Return a unique ID for the intersection point.
-        // Return a positive id for a point, or a negative id for an edge.
-        int id() const { const IntersectionReference &r = ipref(); return (r.point_id >= 0) ? r.point_id : - r.edge_id; }
-        bool operator==(const OpenPolylineEnd &rhs) const { return this->polyline == rhs.polyline && this->start == rhs.start; }
-    };
-    auto by_id_lower = [](const OpenPolylineEnd &ope1, const OpenPolylineEnd &ope2) { return ope1.id() < ope2.id(); };
-    std::vector<OpenPolylineEnd> by_id;
-    by_id.reserve(2 * open_polylines.size());
-    for (OpenPolyline &opl : open_polylines) {
-        if (opl.start.point_id != -1 || opl.start.edge_id != -1)
-            by_id.emplace_back(OpenPolylineEnd(&opl, true));
-        if (try_connect_reversed && (opl.end.point_id != -1 || opl.end.edge_id != -1))
-            by_id.emplace_back(OpenPolylineEnd(&opl, false));
-    }
-    std::sort(by_id.begin(), by_id.end(), by_id_lower);
-    // Find an iterator to by_id_lower for the particular end of OpenPolyline (by comparing the OpenPolyline pointer and the start attribute).
-    auto find_polyline_end = [&by_id, by_id_lower](const OpenPolylineEnd &end) -> std::vector<OpenPolylineEnd>::iterator {
-        for (auto it = std::lower_bound(by_id.begin(), by_id.end(), end, by_id_lower);
-                  it != by_id.end() && it->id() == end.id(); ++ it)
-            if (*it == end)
-                return it;
-        return by_id.end();
-    };
-    // Try to connect the loops.
-    std::vector<OpenPolyline*> sorted_by_length = open_polylines_sorted(open_polylines, false);
-    for (OpenPolyline *opl : sorted_by_length) {
-        if (opl->consumed)
-            continue;
-        opl->consumed = true;
-        OpenPolylineEnd end(opl, false);
-        for (;;) {
-            // find a line starting where last one finishes
-            auto it_next_start = std::lower_bound(by_id.begin(), by_id.end(), end, by_id_lower);
-            for (; it_next_start != by_id.end() && it_next_start->id() == end.id(); ++ it_next_start)
-                if (! it_next_start->polyline->consumed)
-                    goto found;
-            // The current loop could not be closed. Unmark the segment.
-            opl->consumed = false;
-            break;
-        found:
-            // Attach this polyline to the end of the initial polyline.
-            if (it_next_start->start) {
-                auto it = it_next_start->polyline->points.begin();
-                std::copy(++ it, it_next_start->polyline->points.end(), back_inserter(opl->points));
-            } else {
-                auto it = it_next_start->polyline->points.rbegin();
-                std::copy(++ it, it_next_start->polyline->points.rend(), back_inserter(opl->points));
-            }
-            opl->length += it_next_start->polyline->length;
-            // Mark the next polyline as consumed.
-            it_next_start->polyline->points.clear();
-            it_next_start->polyline->length = 0.;
-            it_next_start->polyline->consumed = true;
-            if (try_connect_reversed) {
-                // Running in a mode, where the polylines may be connected by mixing their orientations.
-                // Update the end point lookup structure after the end point of the current polyline was extended.
-                auto it_end      = find_polyline_end(end);
-                auto it_next_end = find_polyline_end(OpenPolylineEnd(it_next_start->polyline, !it_next_start->start));
-                // Swap the end points of the current and next polyline, but keep the polyline ptr and the start flag.
-                std::swap(opl->end, it_next_end->start ? it_next_end->polyline->start : it_next_end->polyline->end);
-                // Swap the positions of OpenPolylineEnd structures in the sorted array to match their respective end point positions.
-                std::swap(*it_end, *it_next_end);
-            }
-            // Check whether we closed this loop.
-            if ((opl->start.edge_id  != -1 && opl->start.edge_id  == opl->end.edge_id) ||
-                (opl->start.point_id != -1 && opl->start.point_id == opl->end.point_id)) {
-                // The current loop is complete. Add it to the output.
-                //assert(opl->points.front().point_id == opl->points.back().point_id);
-                //assert(opl->points.front().edge_id  == opl->points.back().edge_id);
-                // Remove the duplicate last point.
-                opl->points.pop_back();
-                if (opl->points.size() >= 3) {
-                    if (try_connect_reversed && area(opl->points) < 0)
-                        // The closed polygon is patched from pieces with messed up orientation, therefore
-                        // the orientation of the patched up polygon is not known.
-                        // Orient the patched up polygons CCW. This heuristic may close some holes and cavities.
-                        std::reverse(opl->points.begin(), opl->points.end());
-                    loops.emplace_back(std::move(opl->points));
-                }
-                opl->points.clear();
-                break;
-            }
-            // Continue with the current loop.
-        }
-    }
+    return create_face_neighbors_index_impl(its, [](){});
 }
 
-// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices, 
-// possibly closing small gaps.
-// Depending on "try_connect_reversed", it may or may not connect segments crossing triangles of opposite orientation.
-static void chain_open_polylines_close_gaps(std::vector<OpenPolyline> &open_polylines, Polygons &loops, double max_gap, bool try_connect_reversed)
+std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its, std::function<void()> throw_on_cancel_callback)
 {
-    const coord_t max_gap_scaled = (coord_t)scale_(max_gap);
-
-    // Sort the open polylines by their length, so the new loops will be seeded from longer chains.
-    // Update the polyline lengths, return only not yet consumed polylines.
-    std::vector<OpenPolyline*> sorted_by_length = open_polylines_sorted(open_polylines, true);
-
-    // Store the end points of open_polylines into ClosestPointInRadiusLookup<OpenPolylineEnd>.
-    struct OpenPolylineEnd {
-        OpenPolylineEnd(OpenPolyline *polyline, bool start) : polyline(polyline), start(start) {}
-        OpenPolyline    *polyline;
-        // Is it the start or end point?
-        bool             start;
-        const Point&     point() const { return start ? polyline->points.front() : polyline->points.back(); }
-        bool operator==(const OpenPolylineEnd &rhs) const { return this->polyline == rhs.polyline && this->start == rhs.start; }
-    };
-    struct OpenPolylineEndAccessor {
-        const Point* operator()(const OpenPolylineEnd &pt) const { return pt.polyline->consumed ? nullptr : &pt.point(); }
-    };
-    typedef ClosestPointInRadiusLookup<OpenPolylineEnd, OpenPolylineEndAccessor> ClosestPointLookupType;
-    ClosestPointLookupType closest_end_point_lookup(max_gap_scaled);
-    for (OpenPolyline *opl : sorted_by_length) {
-        closest_end_point_lookup.insert(OpenPolylineEnd(opl, true));
-        if (try_connect_reversed)
-            closest_end_point_lookup.insert(OpenPolylineEnd(opl, false));
-    }
-    // Try to connect the loops.
-    for (OpenPolyline *opl : sorted_by_length) {
-        if (opl->consumed)
-            continue;
-        OpenPolylineEnd end(opl, false);
-        if (try_connect_reversed)
-            // The end point of this polyline will be modified, thus the following entry will become invalid. Remove it.
-            closest_end_point_lookup.erase(end);
-        opl->consumed = true;
-        size_t n_segments_joined = 1;
-        for (;;) {
-            // Find a line starting where last one finishes, only return non-consumed open polylines (OpenPolylineEndAccessor returns null for consumed).
-            std::pair<const OpenPolylineEnd*, double> next_start_and_dist = closest_end_point_lookup.find(end.point());
-            const OpenPolylineEnd *next_start = next_start_and_dist.first;
-            // Check whether we closed this loop.
-            double current_loop_closing_distance2 = (opl->points.back() - opl->points.front()).cast<double>().squaredNorm();
-            bool   loop_closed = current_loop_closing_distance2 < coordf_t(max_gap_scaled) * coordf_t(max_gap_scaled);
-            if (next_start != nullptr && loop_closed && current_loop_closing_distance2 < next_start_and_dist.second) {
-                // Heuristics to decide, whether to close the loop, or connect another polyline.
-                // One should avoid closing loops shorter than max_gap_scaled.
-                loop_closed = sqrt(current_loop_closing_distance2) < 0.3 * length(opl->points);
-            }
-            if (loop_closed) {
-                // Remove the start point of the current polyline from the lookup.
-                // Mark the current segment as not consumed, otherwise the closest_end_point_lookup.erase() would fail.
-                opl->consumed = false;
-                closest_end_point_lookup.erase(OpenPolylineEnd(opl, true));
-                if (current_loop_closing_distance2 == 0.) {
-                    // Remove the duplicate last point.
-                    opl->points.pop_back();
-                } else {
-                    // The end points are different, keep both of them.
-                }
-                if (opl->points.size() >= 3) {
-                    if (try_connect_reversed && n_segments_joined > 1 && area(opl->points) < 0)
-                        // The closed polygon is patched from pieces with messed up orientation, therefore
-                        // the orientation of the patched up polygon is not known.
-                        // Orient the patched up polygons CCW. This heuristic may close some holes and cavities.
-                        std::reverse(opl->points.begin(), opl->points.end());
-                    loops.emplace_back(std::move(opl->points));
-                }
-                opl->points.clear();
-                opl->consumed = true;
-                break;
-            }
-            if (next_start == nullptr) {
-                // The current loop could not be closed. Unmark the segment.
-                opl->consumed = false;
-                if (try_connect_reversed)
-                    // Re-insert the end point.
-                    closest_end_point_lookup.insert(OpenPolylineEnd(opl, false));
-                break;
-            }
-            // Attach this polyline to the end of the initial polyline.
-            if (next_start->start) {
-                auto it = next_start->polyline->points.begin();
-                if (*it == opl->points.back())
-                    ++ it;
-                std::copy(it, next_start->polyline->points.end(), back_inserter(opl->points));
-            } else {
-                auto it = next_start->polyline->points.rbegin();
-                if (*it == opl->points.back())
-                    ++ it;
-                std::copy(it, next_start->polyline->points.rend(), back_inserter(opl->points));
-            }
-            ++ n_segments_joined;
-            // Remove the end points of the consumed polyline segment from the lookup.
-            OpenPolyline *opl2 = next_start->polyline;
-            closest_end_point_lookup.erase(OpenPolylineEnd(opl2, true));
-            if (try_connect_reversed)
-                closest_end_point_lookup.erase(OpenPolylineEnd(opl2, false));
-            opl2->points.clear();
-            opl2->consumed = true;
-            // Continue with the current loop.
-        }
-    }
+    return create_face_neighbors_index_impl(its, throw_on_cancel_callback);
 }
 
-void TriangleMeshSlicer::make_loops(std::vector<IntersectionLine> &lines, Polygons* loops) const
+int its_remove_degenerate_faces(indexed_triangle_set &its, bool shrink_to_fit)
 {
-#if 0
-//FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
-//#ifdef _DEBUG
-    for (const Line &l : lines)
-        assert(l.a != l.b);
-#endif /* _DEBUG */
-
-    // There should be no tangent edges, as the horizontal triangles are ignored and if two triangles touch at a cutting plane,
-    // only the bottom triangle is considered to be cutting the plane.
-//    remove_tangent_edges(lines);
-
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-        BoundingBox bbox_svg;
-        {
-            static int iRun = 0;
-            for (const Line &line : lines) {
-                bbox_svg.merge(line.a);
-                bbox_svg.merge(line.b);
-            }
-            SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-raw_lines-%d.svg", iRun ++).c_str(), bbox_svg);
-            for (const Line &line : lines)
-                svg.draw(line);
-            svg.Close();
+    int last = 0;
+    for (int i = 0; i < int(its.indices.size()); ++ i) {
+        const stl_triangle_vertex_indices &face = its.indices[i];
+        if (face(0) != face(1) && face(0) != face(2) && face(1) != face(2)) {
+            if (last < i)
+                its.indices[last] = its.indices[i];
+            ++ last;
         }
-#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-
-    std::vector<OpenPolyline> open_polylines;
-    chain_lines_by_triangle_connectivity(lines, *loops, open_polylines);
-
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-        {
-            static int iRun = 0;
-            SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines-%d.svg", iRun ++).c_str(), bbox_svg);
-            svg.draw(union_ex(*loops));
-            for (const OpenPolyline &pl : open_polylines)
-                svg.draw(Polyline(pl.points), "red");
-            svg.Close();
-        }
-#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-
-    // Now process the open polylines.
-    // Do it in two rounds, first try to connect in the same direction only,
-    // then try to connect the open polylines in reversed order as well.
-    chain_open_polylines_exact(open_polylines, *loops, false);
-    chain_open_polylines_exact(open_polylines, *loops, true);
-
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-    {
-        static int iRun = 0;
-        SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines2-%d.svg", iRun++).c_str(), bbox_svg);
-        svg.draw(union_ex(*loops));
-        for (const OpenPolyline &pl : open_polylines) {
-            if (pl.points.empty())
-                continue;
-            svg.draw(Polyline(pl.points), "red");
-            svg.draw(pl.points.front(), "blue");
-            svg.draw(pl.points.back(), "blue");
-        }
-        svg.Close();
     }
-#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-
-    // Try to close gaps.
-    // Do it in two rounds, first try to connect in the same direction only,
-    // then try to connect the open polylines in reversed order as well.
-#if 0
-    for (double max_gap : { EPSILON, 0.001, 0.1, 1., 2. }) {
-        chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, false);
-        chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, true);
+    int removed = int(its.indices.size()) - last;
+    if (removed) {
+        its.indices.erase(its.indices.begin() + last, its.indices.end());
+        // Optionally shrink the vertices.
+        if (shrink_to_fit)
+            its.indices.shrink_to_fit();
     }
-#else
-    const double max_gap = 2.; //mm
-    chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, false);
-    chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, true);
-#endif
-
-#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-    {
-        static int iRun = 0;
-        SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines-final-%d.svg", iRun++).c_str(), bbox_svg);
-        svg.draw(union_ex(*loops));
-        for (const OpenPolyline &pl : open_polylines) {
-            if (pl.points.empty())
-                continue;
-            svg.draw(Polyline(pl.points), "red");
-            svg.draw(pl.points.front(), "blue");
-            svg.draw(pl.points.back(), "blue");
-        }
-        svg.Close();
-    }
-#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+    return removed;
 }
 
-// Only used to cut the mesh into two halves.
-void TriangleMeshSlicer::make_expolygons_simple(std::vector<IntersectionLine> &lines, ExPolygons* slices) const
+int its_compactify_vertices(indexed_triangle_set &its, bool shrink_to_fit)
 {
-    assert(slices->empty());
-
-    Polygons loops;
-    this->make_loops(lines, &loops);
-    
-    Polygons holes;
-    for (Polygons::const_iterator loop = loops.begin(); loop != loops.end(); ++ loop) {
-        if (loop->area() >= 0.) {
-            ExPolygon ex;
-            ex.contour = *loop;
-            slices->emplace_back(ex);
-        } else {
-            holes.emplace_back(*loop);
+    // First used to mark referenced vertices, later used for mapping old vertex index to a new one.
+    std::vector<int> vertex_map(its.vertices.size(), 0);
+    // Mark referenced vertices.
+    for (const stl_triangle_vertex_indices &face : its.indices)
+        for (int i = 0; i < 3; ++ i)
+            vertex_map[face(i)] = 1;
+    // Compactify vertices, update map from old vertex index to a new one.
+    int last = 0;
+    for (int i = 0; i < int(vertex_map.size()); ++ i)
+        if (vertex_map[i]) {
+            if (last < i)
+                its.vertices[last] = its.vertices[i];
+            vertex_map[i] = last ++;
         }
+    int removed = int(its.vertices.size()) - last;
+    if (removed) {
+        its.vertices.erase(its.vertices.begin() + last, its.vertices.end());
+        // Update faces with the new vertex indices.
+        for (stl_triangle_vertex_indices &face : its.indices)
+            for (int i = 0; i < 3; ++ i)
+                face(i) = vertex_map[face(i)];
+        // Optionally shrink the vertices.
+        if (shrink_to_fit)
+            its.vertices.shrink_to_fit();
     }
-
-    // If there are holes, then there should also be outer contours.
-    assert(holes.empty() || ! slices->empty());
-    if (slices->empty())
-        return;
-    
-    // Assign holes to outer contours.
-    for (Polygons::const_iterator hole = holes.begin(); hole != holes.end(); ++ hole) {
-        // Find an outer contour to a hole.
-        int     slice_idx            = -1;
-        double  current_contour_area = std::numeric_limits<double>::max();
-        for (ExPolygons::iterator slice = slices->begin(); slice != slices->end(); ++ slice) {
-            if (slice->contour.contains(hole->points.front())) {
-                double area = slice->contour.area();
-                if (area < current_contour_area) {
-                    slice_idx = slice - slices->begin();
-                    current_contour_area = area;
-                }
-            }
-        }
-        // assert(slice_idx != -1);
-        if (slice_idx == -1)
-            // Ignore this hole.
-            continue;
-        assert(current_contour_area < std::numeric_limits<double>::max() && current_contour_area >= -hole->area());
-        (*slices)[slice_idx].holes.emplace_back(std::move(*hole));
-    }
-
-#if 0
-    // If the input mesh is not valid, the holes may intersect with the external contour.
-    // Rather subtract them from the outer contour.
-    Polygons poly;
-    for (auto it_slice = slices->begin(); it_slice != slices->end(); ++ it_slice) {
-        if (it_slice->holes.empty()) {
-            poly.emplace_back(std::move(it_slice->contour));
-        } else {
-            Polygons contours;
-            contours.emplace_back(std::move(it_slice->contour));
-            for (auto it = it_slice->holes.begin(); it != it_slice->holes.end(); ++ it)
-                it->reverse();
-            polygons_append(poly, diff(contours, it_slice->holes));
-        }
-    }
-    // If the input mesh is not valid, the input contours may intersect.
-    *slices = union_ex(poly);
-#endif
-
-#if 0
-    // If the input mesh is not valid, the holes may intersect with the external contour.
-    // Rather subtract them from the outer contour.
-    ExPolygons poly;
-    for (auto it_slice = slices->begin(); it_slice != slices->end(); ++ it_slice) {
-        Polygons contours;
-        contours.emplace_back(std::move(it_slice->contour));
-        for (auto it = it_slice->holes.begin(); it != it_slice->holes.end(); ++ it)
-            it->reverse();
-        expolygons_append(poly, diff_ex(contours, it_slice->holes));
-    }
-    // If the input mesh is not valid, the input contours may intersect.
-    *slices = std::move(poly);
-#endif
+    return removed;
 }
 
-void TriangleMeshSlicer::make_expolygons(const Polygons &loops, const float closing_radius, ExPolygons* slices) const
+void its_shrink_to_fit(indexed_triangle_set &its)
 {
-    /*
-        Input loops are not suitable for evenodd nor nonzero fill types, as we might get
-        two consecutive concentric loops having the same winding order - and we have to 
-        respect such order. In that case, evenodd would create wrong inversions, and nonzero
-        would ignore holes inside two concentric contours.
-        So we're ordering loops and collapse consecutive concentric loops having the same 
-        winding order.
-        TODO: find a faster algorithm for this, maybe with some sort of binary search.
-        If we computed a "nesting tree" we could also just remove the consecutive loops
-        having the same winding order, and remove the extra one(s) so that we could just
-        supply everything to offset() instead of performing several union/diff calls.
-    
-        we sort by area assuming that the outermost loops have larger area;
-        the previous sorting method, based on $b->contains($a->[0]), failed to nest
-        loops correctly in some edge cases when original model had overlapping facets
-    */
-
-    /* The following lines are commented out because they can generate wrong polygons,
-       see for example issue #661 */
-
-    //std::vector<double> area;
-    //std::vector<size_t> sorted_area;  // vector of indices
-    //for (Polygons::const_iterator loop = loops.begin(); loop != loops.end(); ++ loop) {
-    //    area.emplace_back(loop->area());
-    //    sorted_area.emplace_back(loop - loops.begin());
-    //}
-    //
-    //// outer first
-    //std::sort(sorted_area.begin(), sorted_area.end(),
-    //    [&area](size_t a, size_t b) { return std::abs(area[a]) > std::abs(area[b]); });
-
-    //// we don't perform a safety offset now because it might reverse cw loops
-    //Polygons p_slices;
-    //for (std::vector<size_t>::const_iterator loop_idx = sorted_area.begin(); loop_idx != sorted_area.end(); ++ loop_idx) {
-    //    /* we rely on the already computed area to determine the winding order
-    //       of the loops, since the Orientation() function provided by Clipper
-    //       would do the same, thus repeating the calculation */
-    //    Polygons::const_iterator loop = loops.begin() + *loop_idx;
-    //    if (area[*loop_idx] > +EPSILON)
-    //        p_slices.emplace_back(*loop);
-    //    else if (area[*loop_idx] < -EPSILON)
-    //        //FIXME This is arbitrary and possibly very slow.
-    //        // If the hole is inside a polygon, then there is no need to diff.
-    //        // If the hole intersects a polygon boundary, then diff it, but then
-    //        // there is no guarantee of an ordering of the loops.
-    //        // Maybe we can test for the intersection before running the expensive diff algorithm?
-    //        p_slices = diff(p_slices, *loop);
-    //}
-
-    // Perform a safety offset to merge very close facets (TODO: find test case for this)
-    // 0.0499 comes from https://github.com/slic3r/Slic3r/issues/959
-//    double safety_offset = scale_(0.0499);
-    // 0.0001 is set to satisfy GH #520, #1029, #1364
-    double safety_offset = scale_(closing_radius);
-
-    /* The following line is commented out because it can generate wrong polygons,
-       see for example issue #661 */
-    //ExPolygons ex_slices = offset2_ex(p_slices, +safety_offset, -safety_offset);
-    
-    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
-    size_t holes_count = 0;
-    for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++ e)
-        holes_count += e->holes.size();
-    printf("%zu surface(s) having %zu holes detected from %zu polylines\n",
-        ex_slices.size(), holes_count, loops.size());
-    #endif
-    
-    // append to the supplied collection
-    if (safety_offset > 0)
-        expolygons_append(*slices, offset2_ex(union_ex(loops), +safety_offset, -safety_offset));
-    else
-        expolygons_append(*slices, union_ex(loops));
-}
-
-void TriangleMeshSlicer::make_expolygons(std::vector<IntersectionLine> &lines, const float closing_radius, ExPolygons* slices) const
-{
-    Polygons pp;
-    this->make_loops(lines, &pp);
-    this->make_expolygons(pp, closing_radius, slices);
-}
-
-void TriangleMeshSlicer::cut(float z, TriangleMesh* upper, TriangleMesh* lower) const
-{
-    IntersectionLines upper_lines, lower_lines;
-    
-    BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - slicing object";
-    float scaled_z = scale_(z);
-    for (uint32_t facet_idx = 0; facet_idx < this->mesh->stl.stats.number_of_facets; ++ facet_idx) {
-        const stl_facet* facet = &this->mesh->stl.facet_start[facet_idx];
-        
-        // find facet extents
-        float min_z = std::min(facet->vertex[0](2), std::min(facet->vertex[1](2), facet->vertex[2](2)));
-        float max_z = std::max(facet->vertex[0](2), std::max(facet->vertex[1](2), facet->vertex[2](2)));
-        
-        // intersect facet with cutting plane
-        IntersectionLine line;
-        if (this->slice_facet(scaled_z, *facet, facet_idx, min_z, max_z, &line) != TriangleMeshSlicer::NoSlice) {
-            // Save intersection lines for generating correct triangulations.
-            if (line.edge_type == feTop) {
-                lower_lines.emplace_back(line);
-            } else if (line.edge_type == feBottom) {
-                upper_lines.emplace_back(line);
-            } else if (line.edge_type != feHorizontal) {
-                lower_lines.emplace_back(line);
-                upper_lines.emplace_back(line);
-            }
-        }
-        
-        if (min_z > z || (min_z == z && max_z > z)) {
-            // facet is above the cut plane and does not belong to it
-            if (upper != nullptr)
-				stl_add_facet(&upper->stl, facet);
-        } else if (max_z < z || (max_z == z && min_z < z)) {
-            // facet is below the cut plane and does not belong to it
-            if (lower != nullptr)
-				stl_add_facet(&lower->stl, facet);
-        } else if (min_z < z && max_z > z) {
-            // Facet is cut by the slicing plane.
-
-            // look for the vertex on whose side of the slicing plane there are no other vertices
-            int isolated_vertex;
-            if ( (facet->vertex[0](2) > z) == (facet->vertex[1](2) > z) ) {
-                isolated_vertex = 2;
-            } else if ( (facet->vertex[1](2) > z) == (facet->vertex[2](2) > z) ) {
-                isolated_vertex = 0;
-            } else {
-                isolated_vertex = 1;
-            }
-            
-            // get vertices starting from the isolated one
-            const stl_vertex &v0 = facet->vertex[isolated_vertex];
-            const stl_vertex &v1 = facet->vertex[(isolated_vertex+1) % 3];
-            const stl_vertex &v2 = facet->vertex[(isolated_vertex+2) % 3];
-            
-            // intersect v0-v1 and v2-v0 with cutting plane and make new vertices
-            stl_vertex v0v1, v2v0;
-            v0v1(0) = v1(0) + (v0(0) - v1(0)) * (z - v1(2)) / (v0(2) - v1(2));
-            v0v1(1) = v1(1) + (v0(1) - v1(1)) * (z - v1(2)) / (v0(2) - v1(2));
-            v0v1(2) = z;
-            v2v0(0) = v2(0) + (v0(0) - v2(0)) * (z - v2(2)) / (v0(2) - v2(2));
-            v2v0(1) = v2(1) + (v0(1) - v2(1)) * (z - v2(2)) / (v0(2) - v2(2));
-            v2v0(2) = z;
-            
-            // build the triangular facet
-            stl_facet triangle;
-            triangle.normal = facet->normal;
-            triangle.vertex[0] = v0;
-            triangle.vertex[1] = v0v1;
-            triangle.vertex[2] = v2v0;
-            
-            // build the facets forming a quadrilateral on the other side
-            stl_facet quadrilateral[2];
-            quadrilateral[0].normal = facet->normal;
-            quadrilateral[0].vertex[0] = v1;
-            quadrilateral[0].vertex[1] = v2;
-            quadrilateral[0].vertex[2] = v0v1;
-            quadrilateral[1].normal = facet->normal;
-            quadrilateral[1].vertex[0] = v2;
-            quadrilateral[1].vertex[1] = v2v0;
-            quadrilateral[1].vertex[2] = v0v1;
-            
-            if (v0(2) > z) {
-                if (upper != nullptr) 
-					stl_add_facet(&upper->stl, &triangle);
-                if (lower != nullptr) {
-                    stl_add_facet(&lower->stl, &quadrilateral[0]);
-                    stl_add_facet(&lower->stl, &quadrilateral[1]);
-                }
-            } else {
-                if (upper != nullptr) {
-                    stl_add_facet(&upper->stl, &quadrilateral[0]);
-                    stl_add_facet(&upper->stl, &quadrilateral[1]);
-                }
-                if (lower != nullptr) 
-					stl_add_facet(&lower->stl, &triangle);
-            }
-        }
-    }
-    
-    if (upper != nullptr) {
-        BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - triangulating upper part";
-        ExPolygons section;
-        this->make_expolygons_simple(upper_lines, &section);
-        Pointf3s triangles = triangulate_expolygons_3d(section, z, true);
-        stl_facet facet;
-        facet.normal = stl_normal(0, 0, -1.f);
-        for (size_t i = 0; i < triangles.size(); ) {
-            for (size_t j = 0; j < 3; ++ j)
-                facet.vertex[j] = triangles[i ++].cast<float>();
-            stl_add_facet(&upper->stl, &facet);
-        }
-    }
-    
-    if (lower != nullptr) {
-        BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - triangulating lower part";
-        ExPolygons section;
-        this->make_expolygons_simple(lower_lines, &section);
-        Pointf3s triangles = triangulate_expolygons_3d(section, z, false);
-        stl_facet facet;
-        facet.normal = stl_normal(0, 0, -1.f);
-        for (size_t i = 0; i < triangles.size(); ) {
-            for (size_t j = 0; j < 3; ++ j)
-                facet.vertex[j] = triangles[i ++].cast<float>();
-            stl_add_facet(&lower->stl, &facet);
-        }
-    }
-    
-    BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - updating object sizes";
-    stl_get_size(&upper->stl);
-    stl_get_size(&lower->stl);
+    its.indices.shrink_to_fit();
+    its.vertices.shrink_to_fit();
 }
 
 // Generate the vertex list for a cube solid of arbitrary size in X/Y/Z.
@@ -1962,8 +851,8 @@ TriangleMesh make_cube(double x, double y, double z)
             {1, 7, 6}, {1, 6, 2}, {2, 6, 5},
             {2, 5, 3}, {4, 0, 3}, {4, 3, 5}
         });
-	mesh.repair();
-	return mesh;
+    mesh.repair();
+    return mesh;
 }
 
 // Generate the mesh for a cylinder and return it, using 
@@ -1971,13 +860,13 @@ TriangleMesh make_cube(double x, double y, double z)
 // Default is 360 sides, angle fa is in radians.
 TriangleMesh make_cylinder(double r, double h, double fa)
 {
-	size_t n_steps    = (size_t)ceil(2. * PI / fa);
-	double angle_step = 2. * PI / n_steps;
+    size_t n_steps    = (size_t)ceil(2. * PI / fa);
+    double angle_step = 2. * PI / n_steps;
 
-	Pointf3s			vertices;
-	std::vector<Vec3i>	facets;
-	vertices.reserve(2 * n_steps + 2);
-	facets.reserve(4 * n_steps);
+    Pointf3s            vertices;
+    std::vector<Vec3i>  facets;
+    vertices.reserve(2 * n_steps + 2);
+    facets.reserve(4 * n_steps);
 
     // 2 special vertices, top and bottom center, rest are relative to this
     vertices.emplace_back(Vec3d(0.0, 0.0, 0.0));
@@ -1987,29 +876,29 @@ TriangleMesh make_cylinder(double r, double h, double fa)
     // circle, generate four points and four facets (2 for the wall, 2 for the
     // top and bottom.
     // Special case: Last line shares 2 vertices with the first line.
-	Vec2d p = Eigen::Rotation2Dd(0.) * Eigen::Vector2d(0, r);
-	vertices.emplace_back(Vec3d(p(0), p(1), 0.));
-	vertices.emplace_back(Vec3d(p(0), p(1), h));
-	for (size_t i = 1; i < n_steps; ++i) {
+    Vec2d p = Eigen::Rotation2Dd(0.) * Eigen::Vector2d(0, r);
+    vertices.emplace_back(Vec3d(p(0), p(1), 0.));
+    vertices.emplace_back(Vec3d(p(0), p(1), h));
+    for (size_t i = 1; i < n_steps; ++i) {
         p = Eigen::Rotation2Dd(angle_step * i) * Eigen::Vector2d(0, r);
         vertices.emplace_back(Vec3d(p(0), p(1), 0.));
         vertices.emplace_back(Vec3d(p(0), p(1), h));
         int id = (int)vertices.size() - 1;
         facets.emplace_back( 0, id - 1, id - 3); // top
         facets.emplace_back(id,      1, id - 2); // bottom
-		facets.emplace_back(id, id - 2, id - 3); // upper-right of side
+        facets.emplace_back(id, id - 2, id - 3); // upper-right of side
         facets.emplace_back(id, id - 3, id - 1); // bottom-left of side
     }
     // Connect the last set of vertices with the first.
-	int id = (int)vertices.size() - 1;
+    int id = (int)vertices.size() - 1;
     facets.emplace_back( 0, 2, id - 1);
     facets.emplace_back( 3, 1,     id);
-	facets.emplace_back(id, 2,      3);
+    facets.emplace_back(id, 2,      3);
     facets.emplace_back(id, id - 1, 2);
     
-	TriangleMesh mesh(std::move(vertices), std::move(facets));
-	mesh.repair();
-	return mesh;
+    TriangleMesh mesh(std::move(vertices), std::move(facets));
+    mesh.repair();
+    return mesh;
 }
 
 // Generates mesh for a sphere centered about the origin, using the generated angle
@@ -2018,74 +907,56 @@ TriangleMesh make_cylinder(double r, double h, double fa)
 //FIXME better to discretize an Icosahedron recursively http://www.songho.ca/opengl/gl_sphere.html
 TriangleMesh make_sphere(double radius, double fa)
 {
-	int   sectorCount = int(ceil(2. * M_PI / fa));
-	int   stackCount  = int(ceil(M_PI / fa));
-	float sectorStep  = float(2. * M_PI / sectorCount);
-	float stackStep   = float(M_PI / stackCount);
+    int   sectorCount = int(ceil(2. * M_PI / fa));
+    int   stackCount  = int(ceil(M_PI / fa));
+    float sectorStep  = float(2. * M_PI / sectorCount);
+    float stackStep   = float(M_PI / stackCount);
 
-	Pointf3s vertices;
-	vertices.reserve((stackCount - 1) * sectorCount + 2);
-	for (int i = 0; i <= stackCount; ++ i) {
-		// from pi/2 to -pi/2
-		double stackAngle = 0.5 * M_PI - stackStep * i;
-		double xy = radius * cos(stackAngle);
-		double z  = radius * sin(stackAngle);
-		if (i == 0 || i == stackCount)
-			vertices.emplace_back(Vec3d(xy, 0., z));
-		else
-			for (int j = 0; j < sectorCount; ++ j) {
-				// from 0 to 2pi
-				double sectorAngle = sectorStep * j;
-				vertices.emplace_back(Vec3d(xy * cos(sectorAngle), xy * sin(sectorAngle), z));
-			}
-	}
-
-	std::vector<Vec3i> facets;
-	facets.reserve(2 * (stackCount - 1) * sectorCount);
-	for (int i = 0; i < stackCount; ++ i) {
-		// Beginning of current stack.
-		int k1 = (i == 0) ? 0 : (1 + (i - 1) * sectorCount);
-		int k1_first = k1;
-		// Beginning of next stack.
-		int k2 = (i == 0) ? 1 : (k1 + sectorCount);
-		int k2_first = k2;
-		for (int j = 0; j < sectorCount; ++ j) {
-			// 2 triangles per sector excluding first and last stacks
-			int k1_next = k1;
-			int k2_next = k2;
-			if (i != 0) {
-				k1_next = (j + 1 == sectorCount) ? k1_first : (k1 + 1);
-				facets.emplace_back(k1, k2, k1_next);
-			}
-			if (i + 1 != stackCount) {
-				k2_next = (j + 1 == sectorCount) ? k2_first : (k2 + 1);
-				facets.emplace_back(k1_next, k2, k2_next);
-			}
-			k1 = k1_next;
-			k2 = k2_next;
-		}
-	}
-	TriangleMesh mesh(std::move(vertices), std::move(facets));
-	mesh.repair();
-	return mesh;
-}
-
-std::vector<std::vector<size_t> > create_neighbor_index(const indexed_triangle_set &its)
-{
-    if (its.vertices.empty()) return {};
-
-    size_t res = its.indices.size() / its.vertices.size();
-    std::vector< std::vector<size_t> > index(its.vertices.size(),
-                                             reserve_vector<size_t>(res));
-
-    for (size_t fi = 0; fi < its.indices.size(); ++fi) {
-        auto &face = its.indices[fi];
-        index[face(0)].emplace_back(fi);
-        index[face(1)].emplace_back(fi);
-        index[face(2)].emplace_back(fi);
+    Pointf3s vertices;
+    vertices.reserve((stackCount - 1) * sectorCount + 2);
+    for (int i = 0; i <= stackCount; ++ i) {
+        // from pi/2 to -pi/2
+        double stackAngle = 0.5 * M_PI - stackStep * i;
+        double xy = radius * cos(stackAngle);
+        double z  = radius * sin(stackAngle);
+        if (i == 0 || i == stackCount)
+            vertices.emplace_back(Vec3d(xy, 0., z));
+        else
+            for (int j = 0; j < sectorCount; ++ j) {
+                // from 0 to 2pi
+                double sectorAngle = sectorStep * j;
+                vertices.emplace_back(Vec3d(xy * cos(sectorAngle), xy * sin(sectorAngle), z));
+            }
     }
 
-    return index;
+    std::vector<Vec3i> facets;
+    facets.reserve(2 * (stackCount - 1) * sectorCount);
+    for (int i = 0; i < stackCount; ++ i) {
+        // Beginning of current stack.
+        int k1 = (i == 0) ? 0 : (1 + (i - 1) * sectorCount);
+        int k1_first = k1;
+        // Beginning of next stack.
+        int k2 = (i == 0) ? 1 : (k1 + sectorCount);
+        int k2_first = k2;
+        for (int j = 0; j < sectorCount; ++ j) {
+            // 2 triangles per sector excluding first and last stacks
+            int k1_next = k1;
+            int k2_next = k2;
+            if (i != 0) {
+                k1_next = (j + 1 == sectorCount) ? k1_first : (k1 + 1);
+                facets.emplace_back(k1, k2, k1_next);
+            }
+            if (i + 1 != stackCount) {
+                k2_next = (j + 1 == sectorCount) ? k2_first : (k2 + 1);
+                facets.emplace_back(k1_next, k2, k2_next);
+            }
+            k1 = k1_next;
+            k2 = k2_next;
+        }
+    }
+    TriangleMesh mesh(std::move(vertices), std::move(facets));
+    mesh.repair();
+    return mesh;
 }
 
 }
diff --git a/src/libslic3r/TriangleMesh.hpp b/src/libslic3r/TriangleMesh.hpp
index e6f6dc84b..dd1a76156 100644
--- a/src/libslic3r/TriangleMesh.hpp
+++ b/src/libslic3r/TriangleMesh.hpp
@@ -5,7 +5,6 @@
 #include <admesh/stl.h>
 #include <functional>
 #include <vector>
-#include <boost/thread.hpp>
 #include "BoundingBox.hpp"
 #include "Line.hpp"
 #include "Point.hpp"
@@ -24,7 +23,7 @@ public:
     TriangleMesh() : repaired(false) {}
     TriangleMesh(const Pointf3s &points, const std::vector<Vec3i> &facets);
     explicit TriangleMesh(const indexed_triangle_set &M);
-	void clear() { this->stl.clear(); this->its.clear(); this->repaired = false; }
+    void clear() { this->stl.clear(); this->its.clear(); this->repaired = false; }
     bool ReadSTLFile(const char* input_file) { return stl_open(&stl, input_file); }
     bool write_ascii(const char* output_file) { return stl_write_ascii(&this->stl, output_file, ""); }
     bool write_binary(const char* output_file) { return stl_write_binary(&this->stl, output_file, ""); }
@@ -47,7 +46,7 @@ public:
     void mirror_y() { this->mirror(Y); }
     void mirror_z() { this->mirror(Z); }
     void transform(const Transform3d& t, bool fix_left_handed = false);
-	void transform(const Matrix3d& t, bool fix_left_handed = false);
+    void transform(const Matrix3d& t, bool fix_left_handed = false);
     void align_to_origin();
     void rotate(double angle, Point* center);
     TriangleMeshPtrs split() const;
@@ -62,7 +61,7 @@ public:
     // Return the size of the mesh in coordinates.
     Vec3d size() const { return stl.stats.size.cast<double>(); }
     /// Return the center of the related bounding box.
-	Vec3d center() const { return this->bounding_box().center(); }
+    Vec3d center() const { return this->bounding_box().center(); }
     // Returns the convex hull of this TriangleMesh
     TriangleMesh convex_hull_3d() const;
     // Slice this mesh at the provided Z levels and return the vector
@@ -78,8 +77,8 @@ public:
     size_t memsize() const;
     // Release optional data from the mesh if the object is on the Undo / Redo stack only. Returns the amount of memory released.
     size_t release_optional();
-	// Restore optional data possibly released by release_optional().
-	void restore_optional();
+    // Restore optional data possibly released by release_optional().
+    void restore_optional();
 
     stl_file stl;
     indexed_triangle_set its;
@@ -92,160 +91,16 @@ private:
 // Create an index of faces belonging to each vertex. The returned vector can
 // be indexed with vertex indices and contains a list of face indices for each
 // vertex.
-std::vector< std::vector<size_t> >
-create_neighbor_index(const indexed_triangle_set &its);
+std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangle_set &its);
 
-enum FacetEdgeType { 
-    // A general case, the cutting plane intersect a face at two different edges.
-    feGeneral,
-    // Two vertices are aligned with the cutting plane, the third vertex is below the cutting plane.
-    feTop,
-    // Two vertices are aligned with the cutting plane, the third vertex is above the cutting plane.
-    feBottom,
-    // All three vertices of a face are aligned with the cutting plane.
-    feHorizontal
-};
-
-class IntersectionReference
-{
-public:
-    IntersectionReference() : point_id(-1), edge_id(-1) {}
-    IntersectionReference(int point_id, int edge_id) : point_id(point_id), edge_id(edge_id) {}
-    // Where is this intersection point located? On mesh vertex or mesh edge?
-    // Only one of the following will be set, the other will remain set to -1.
-    // Index of the mesh vertex.
-    int point_id;
-    // Index of the mesh edge.
-    int edge_id;
-};
-
-class IntersectionPoint : public Point, public IntersectionReference
-{
-public:
-    IntersectionPoint() {}
-    IntersectionPoint(int point_id, int edge_id, const Point &pt) : IntersectionReference(point_id, edge_id), Point(pt) {}
-    IntersectionPoint(const IntersectionReference &ir, const Point &pt) : IntersectionReference(ir), Point(pt) {}
-    // Inherits coord_t x, y
-};
-
-class IntersectionLine : public Line
-{
-public:
-    IntersectionLine() : a_id(-1), b_id(-1), edge_a_id(-1), edge_b_id(-1), edge_type(feGeneral), flags(0) {}
-
-    bool skip() const { return (this->flags & SKIP) != 0; }
-    void set_skip() { this->flags |= SKIP; }
-
-    bool is_seed_candidate() const { return (this->flags & NO_SEED) == 0 && ! this->skip(); }
-    void set_no_seed(bool set) { if (set) this->flags |= NO_SEED; else this->flags &= ~NO_SEED; }
-    
-    // Inherits Point a, b
-    // For each line end point, either {a,b}_id or {a,b}edge_a_id is set, the other is left to -1.
-    // Vertex indices of the line end points.
-    int             a_id;
-    int             b_id;
-    // Source mesh edges of the line end points.
-    int             edge_a_id;
-    int             edge_b_id;
-    // feGeneral, feTop, feBottom, feHorizontal
-    FacetEdgeType   edge_type;
-    // Used by TriangleMeshSlicer::slice() to skip duplicate edges.
-    enum {
-        // Triangle edge added, because it has no neighbor.
-        EDGE0_NO_NEIGHBOR   = 0x001,
-        EDGE1_NO_NEIGHBOR   = 0x002,
-        EDGE2_NO_NEIGHBOR   = 0x004,
-        // Triangle edge added, because it makes a fold with another horizontal edge.
-        EDGE0_FOLD          = 0x010,
-        EDGE1_FOLD          = 0x020,
-        EDGE2_FOLD          = 0x040,
-        // The edge cannot be a seed of a greedy loop extraction (folds are not safe to become seeds).
-        NO_SEED             = 0x100,
-        SKIP                = 0x200,
-    };
-    uint32_t        flags;
-};
-typedef std::vector<IntersectionLine> IntersectionLines;
-typedef std::vector<IntersectionLine*> IntersectionLinePtrs;
-
-enum class SlicingMode : uint32_t {
-	// Regular slicing, maintain all contours and their orientation.
-	Regular,
-	// Maintain all contours, orient all contours CCW, therefore all holes are being closed.
-	Positive,
-	// Orient all contours CCW and keep only the contour with the largest area.
-	// This mode is useful for slicing complex objects in vase mode.
-	PositiveLargestContour,
-};
-
-class TriangleMeshSlicer
-{
-public:
-    typedef std::function<void()> throw_on_cancel_callback_type;
-    TriangleMeshSlicer() : mesh(nullptr) {}
-	TriangleMeshSlicer(const TriangleMesh* mesh) { this->init(mesh, [](){}); }
-    void init(const TriangleMesh *mesh, throw_on_cancel_callback_type throw_on_cancel);
-    void slice(
-        const std::vector<float> &z, SlicingMode mode, size_t alternate_mode_first_n_layers, SlicingMode alternate_mode,
-        std::vector<Polygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const;
-    void slice(const std::vector<float> &z, SlicingMode mode, std::vector<Polygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const
-        { return this->slice(z, mode, 0, mode, layers, throw_on_cancel); }
-    void slice(
-        const std::vector<float> &z, SlicingMode mode, size_t alternate_mode_first_n_layers, SlicingMode alternate_mode, const float closing_radius,
-        std::vector<ExPolygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const;
-    void slice(const std::vector<float> &z, SlicingMode mode, const float closing_radius, 
-        std::vector<ExPolygons>* layers, throw_on_cancel_callback_type throw_on_cancel) const
-        { this->slice(z, mode, 0, mode, closing_radius, layers, throw_on_cancel); }
-    enum FacetSliceType {
-        NoSlice = 0,
-        Slicing = 1,
-        Cutting = 2
-    };
-    FacetSliceType slice_facet(float slice_z, const stl_facet &facet, const int facet_idx,
-        const float min_z, const float max_z, IntersectionLine *line_out) const;
-    void cut(float z, TriangleMesh* upper, TriangleMesh* lower) const;
-    void set_up_direction(const Vec3f& up);
-    
-private:
-    const TriangleMesh      *mesh;
-    // Map from a facet to an edge index.
-    std::vector<int>         facets_edges;
-    // Scaled copy of this->mesh->stl.v_shared
-    std::vector<stl_vertex>  v_scaled_shared;
-    // Quaternion that will be used to rotate every facet before the slicing
-    Eigen::Quaternion<float, Eigen::DontAlign> m_quaternion;
-    // Whether or not the above quaterion should be used
-    bool                     m_use_quaternion = false;
-
-    void _slice_do(size_t facet_idx, std::vector<IntersectionLines>* lines, boost::mutex* lines_mutex, const std::vector<float> &z) const;
-    void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops) const;
-    void make_expolygons(const Polygons &loops, const float closing_radius, ExPolygons* slices) const;
-    void make_expolygons_simple(std::vector<IntersectionLine> &lines, ExPolygons* slices) const;
-    void make_expolygons(std::vector<IntersectionLine> &lines, const float closing_radius, ExPolygons* slices) const;
-};
-
-inline void slice_mesh(
-    const TriangleMesh &                              mesh,
-    const std::vector<float> &                        z,
-    std::vector<Polygons> &                           layers,
-    TriangleMeshSlicer::throw_on_cancel_callback_type thr = nullptr)
-{
-    if (mesh.empty()) return;
-    TriangleMeshSlicer slicer(&mesh);
-    slicer.slice(z, SlicingMode::Regular, &layers, thr);
-}
-
-inline void slice_mesh(
-    const TriangleMesh &                              mesh,
-    const std::vector<float> &                        z,
-    std::vector<ExPolygons> &                         layers,
-    float                                             closing_radius,
-    TriangleMeshSlicer::throw_on_cancel_callback_type thr = nullptr)
-{
-    if (mesh.empty()) return;
-    TriangleMeshSlicer slicer(&mesh);
-    slicer.slice(z, SlicingMode::Regular, closing_radius, &layers, thr);
-}
+// Map from a facet edge to a neighbor face index or -1 if no neighbor exists.
+std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its);
+std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its, std::function<void()> throw_on_cancel_callback);
+// Remove degenerate faces, return number of faces removed.
+int its_remove_degenerate_faces(indexed_triangle_set &its, bool shrink_to_fit = true);
+// Remove vertices, which none of the faces references. Return number of freed vertices.
+int its_compactify_vertices(indexed_triangle_set &its, bool shrink_to_fit = true);
+void its_shrink_to_fit(indexed_triangle_set &its);
 
 TriangleMesh make_cube(double x, double y, double z);
 
@@ -259,21 +114,21 @@ TriangleMesh make_sphere(double rho, double fa=(2*PI/360));
 // Serialization through the Cereal library
 #include <cereal/access.hpp>
 namespace cereal {
-	template <class Archive> struct specialize<Archive, Slic3r::TriangleMesh, cereal::specialization::non_member_load_save> {};
-	template<class Archive> void load(Archive &archive, Slic3r::TriangleMesh &mesh) {
+    template <class Archive> struct specialize<Archive, Slic3r::TriangleMesh, cereal::specialization::non_member_load_save> {};
+    template<class Archive> void load(Archive &archive, Slic3r::TriangleMesh &mesh) {
         stl_file &stl = mesh.stl;
         stl.stats.type = inmemory;
-		archive(stl.stats.number_of_facets, stl.stats.original_num_facets);
+        archive(stl.stats.number_of_facets, stl.stats.original_num_facets);
         stl_allocate(&stl);
-		archive.loadBinary((char*)stl.facet_start.data(), stl.facet_start.size() * 50);
+        archive.loadBinary((char*)stl.facet_start.data(), stl.facet_start.size() * 50);
         stl_get_size(&stl);
         mesh.repair();
-	}
-	template<class Archive> void save(Archive &archive, const Slic3r::TriangleMesh &mesh) {
-		const stl_file& stl = mesh.stl;
-		archive(stl.stats.number_of_facets, stl.stats.original_num_facets);
-		archive.saveBinary((char*)stl.facet_start.data(), stl.facet_start.size() * 50);
-	}
+    }
+    template<class Archive> void save(Archive &archive, const Slic3r::TriangleMesh &mesh) {
+        const stl_file& stl = mesh.stl;
+        archive(stl.stats.number_of_facets, stl.stats.original_num_facets);
+        archive.saveBinary((char*)stl.facet_start.data(), stl.facet_start.size() * 50);
+    }
 }
 
 #endif
diff --git a/src/libslic3r/TriangleMeshSlicer.cpp b/src/libslic3r/TriangleMeshSlicer.cpp
new file mode 100644
index 000000000..657a59b62
--- /dev/null
+++ b/src/libslic3r/TriangleMeshSlicer.cpp
@@ -0,0 +1,1363 @@
+#include "ClipperUtils.hpp"
+#include "Geometry.hpp"
+#include "Tesselate.hpp"
+#include "TriangleMesh.hpp"
+#include "TriangleMeshSlicer.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <deque>
+#include <queue>
+#include <utility>
+
+#include <boost/log/trivial.hpp>
+
+#include <tbb/parallel_for.h>
+
+#if 0
+    #define DEBUG
+    #define _DEBUG
+    #undef NDEBUG
+    #define SLIC3R_DEBUG
+// #define SLIC3R_TRIANGLEMESH_DEBUG
+#endif
+
+#include <assert.h>
+
+#if defined(SLIC3R_DEBUG) || defined(SLIC3R_DEBUG_SLICE_PROCESSING)
+#include "SVG.hpp"
+#endif
+
+namespace Slic3r {
+
+class IntersectionReference
+{
+public:
+    IntersectionReference() = default;
+    IntersectionReference(int point_id, int edge_id) : point_id(point_id), edge_id(edge_id) {}
+    // Where is this intersection point located? On mesh vertex or mesh edge?
+    // Only one of the following will be set, the other will remain set to -1.
+    // Index of the mesh vertex.
+    int point_id { -1 };
+    // Index of the mesh edge.
+    int edge_id { -1 };
+};
+
+class IntersectionPoint : public Point, public IntersectionReference
+{
+public:
+    IntersectionPoint() = default;
+    IntersectionPoint(int point_id, int edge_id, const Point &pt) : IntersectionReference(point_id, edge_id), Point(pt) {}
+    IntersectionPoint(const IntersectionReference &ir, const Point &pt) : IntersectionReference(ir), Point(pt) {}
+    // Inherits coord_t x, y
+};
+
+class IntersectionLine : public Line
+{
+public:
+    IntersectionLine() = default;
+
+    bool skip() const { return (this->flags & SKIP) != 0; }
+    void set_skip() { this->flags |= SKIP; }
+
+    bool is_seed_candidate() const { return (this->flags & NO_SEED) == 0 && ! this->skip(); }
+    void set_no_seed(bool set) { if (set) this->flags |= NO_SEED; else this->flags &= ~NO_SEED; }
+    
+    // Inherits Point a, b
+    // For each line end point, either {a,b}_id or {a,b}edge_a_id is set, the other is left to -1.
+    // Vertex indices of the line end points.
+    int             a_id { -1 };
+    int             b_id { -1 };
+    // Source mesh edges of the line end points.
+    int             edge_a_id { -1 };
+    int             edge_b_id { -1 };
+
+    enum class FacetEdgeType { 
+        // A general case, the cutting plane intersect a face at two different edges.
+        General,
+        // Two vertices are aligned with the cutting plane, the third vertex is below the cutting plane.
+        Top,
+        // Two vertices are aligned with the cutting plane, the third vertex is above the cutting plane.
+        Bottom,
+        // All three vertices of a face are aligned with the cutting plane.
+        Horizontal
+    };
+
+    // feGeneral, feTop, feBottom, feHorizontal
+    FacetEdgeType   edge_type { FacetEdgeType::General };
+    // Used by TriangleMeshSlicer::slice() to skip duplicate edges.
+    enum {
+        // Triangle edge added, because it has no neighbor.
+        EDGE0_NO_NEIGHBOR   = 0x001,
+        EDGE1_NO_NEIGHBOR   = 0x002,
+        EDGE2_NO_NEIGHBOR   = 0x004,
+        // Triangle edge added, because it makes a fold with another horizontal edge.
+        EDGE0_FOLD          = 0x010,
+        EDGE1_FOLD          = 0x020,
+        EDGE2_FOLD          = 0x040,
+        // The edge cannot be a seed of a greedy loop extraction (folds are not safe to become seeds).
+        NO_SEED             = 0x100,
+        SKIP                = 0x200,
+    };
+    uint32_t        flags { 0 };
+};
+
+using IntersectionLines = std::vector<IntersectionLine>;
+
+enum class FacetSliceType {
+    NoSlice = 0,
+    Slicing = 1,
+    Cutting = 2
+};
+
+// Return true, if the facet has been sliced and line_out has been filled.
+static FacetSliceType slice_facet(
+    float slice_z, const stl_vertex *vertices, const stl_triangle_vertex_indices &indices, const int *edge_neighbor,
+    const int idx_vertex_lowest, const bool horizontal, IntersectionLine *line_out)
+{
+    IntersectionPoint points[3];
+    size_t            num_points = 0;
+    auto              point_on_layer = size_t(-1);
+
+    // Reorder vertices so that the first one is the one with lowest Z.
+    // This is needed to get all intersection lines in a consistent order
+    // (external on the right of the line)
+    for (int j = 0; j < 3; ++ j) {  // loop through facet edges
+        int               edge_id;
+        const stl_vertex *a, *b;
+        int               a_id, b_id;
+        {
+            int   k = (idx_vertex_lowest + j) % 3;
+            int   l = (k + 1) % 3;
+            edge_id = edge_neighbor[k];
+            a_id    = indices[k];
+            a       = vertices + k;
+            b_id    = indices[l];
+            b       = vertices + l;
+        }
+
+        // Is edge or face aligned with the cutting plane?
+        if (a->z() == slice_z && b->z() == slice_z) {
+            // Edge is horizontal and belongs to the current layer.
+            // The following rotation of the three vertices may not be efficient, but this branch happens rarely.
+            const stl_vertex &v0 = vertices[0];
+            const stl_vertex &v1 = vertices[1];
+            const stl_vertex &v2 = vertices[2];
+            // We may ignore this edge for slicing purposes, but we may still use it for object cutting.
+            FacetSliceType    result = FacetSliceType::Slicing;
+            if (horizontal) {
+                // All three vertices are aligned with slice_z.
+                line_out->edge_type = IntersectionLine::FacetEdgeType::Horizontal;
+                result = FacetSliceType::Cutting;
+                double normal = (v1.x() - v0.x()) * (v2.y() - v1.y()) - (v1.y() - v0.y()) * (v2.x() - v1.x());
+                if (normal < 0) {
+                    // If normal points downwards this is a bottom horizontal facet so we reverse its point order.
+                    std::swap(a, b);
+                    std::swap(a_id, b_id);
+                }
+            } else {
+                // Two vertices are aligned with the cutting plane, the third vertex is below or above the cutting plane.
+                // Is the third vertex below the cutting plane?
+                bool third_below = v0.z() < slice_z || v1.z() < slice_z || v2.z() < slice_z;
+                // Two vertices on the cutting plane, the third vertex is below the plane. Consider the edge to be part of the slice
+                // only if it is the upper edge.
+                // (the bottom most edge resp. vertex of a triangle is not owned by the triangle, but the top most edge resp. vertex is part of the triangle
+                // in respect to the cutting plane).
+                result = third_below ? FacetSliceType::Slicing : FacetSliceType::Cutting;
+                if (third_below) {
+                    line_out->edge_type = IntersectionLine::FacetEdgeType::Top;
+                    std::swap(a, b);
+                    std::swap(a_id, b_id);
+                } else
+                    line_out->edge_type = IntersectionLine::FacetEdgeType::Bottom;
+            }
+            line_out->a.x()  = a->x();
+            line_out->a.y()  = a->y();
+            line_out->b.x()  = b->x();
+            line_out->b.y()  = b->y();
+            line_out->a_id   = a_id;
+            line_out->b_id   = b_id;
+            assert(line_out->a != line_out->b);
+            return result;
+        }
+
+        if (a->z() == slice_z) {
+            // Only point a alings with the cutting plane.
+            if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != a_id) {
+                point_on_layer = num_points;
+                IntersectionPoint &point = points[num_points ++];
+                point.x()      = a->x();
+                point.y()      = a->y();
+                point.point_id = a_id;
+            }
+        } else if (b->z() == slice_z) {
+            // Only point b alings with the cutting plane.
+            if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != b_id) {
+                point_on_layer = num_points;
+                IntersectionPoint &point = points[num_points ++];
+                point.x()      = b->x();
+                point.y()      = b->y();
+                point.point_id = b_id;
+            }
+        } else if ((a->z() < slice_z && b->z() > slice_z) || (b->z() < slice_z && a->z() > slice_z)) {
+            // A general case. The face edge intersects the cutting plane. Calculate the intersection point.
+            assert(a_id != b_id);
+            // Sort the edge to give a consistent answer.
+            if (a_id > b_id) {
+                std::swap(a_id, b_id);
+                std::swap(a, b);
+            }
+            IntersectionPoint &point = points[num_points];
+            double t = (double(slice_z) - double(b->z())) / (double(a->z()) - double(b->z()));
+            if (t <= 0.) {
+                if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != a_id) {
+                    point.x() = a->x();
+                    point.y() = a->y();
+                    point_on_layer = num_points ++;
+                    point.point_id = a_id;
+                }
+            } else if (t >= 1.) {
+                if (point_on_layer == size_t(-1) || points[point_on_layer].point_id != b_id) {
+                    point.x() = b->x();
+                    point.y() = b->y();
+                    point_on_layer = num_points ++;
+                    point.point_id = b_id;
+                }
+            } else {
+                point.x() = coord_t(floor(double(b->x()) + (double(a->x()) - double(b->x())) * t + 0.5));
+                point.y() = coord_t(floor(double(b->y()) + (double(a->y()) - double(b->y())) * t + 0.5));
+                point.edge_id = edge_id;
+                ++ num_points;
+            }
+        }
+    }
+
+    // Facets must intersect each plane 0 or 2 times, or it may touch the plane at a single vertex only.
+    assert(num_points < 3);
+    if (num_points == 2) {
+        line_out->edge_type  = IntersectionLine::FacetEdgeType::General;
+        line_out->a          = (Point)points[1];
+        line_out->b          = (Point)points[0];
+        line_out->a_id       = points[1].point_id;
+        line_out->b_id       = points[0].point_id;
+        line_out->edge_a_id  = points[1].edge_id;
+        line_out->edge_b_id  = points[0].edge_id;
+        // Not a zero lenght edge.
+        //FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
+        //assert(line_out->a != line_out->b);
+        // The plane cuts at least one edge in a general position.
+        assert(line_out->a_id == -1 || line_out->b_id == -1);
+        assert(line_out->edge_a_id != -1 || line_out->edge_b_id != -1);
+        // General slicing position, use the segment for both slicing and object cutting.
+#if 0
+        if (line_out->a_id != -1 && line_out->b_id != -1) {
+            // Solving a degenerate case, where both the intersections snapped to an edge.
+            // Correctly classify the face as below or above based on the position of the 3rd point.
+            int i = indices[0];
+            if (i == line_out->a_id || i == line_out->b_id)
+                i = indices[1];
+            if (i == line_out->a_id || i == line_out->b_id)
+                i = indices[2];
+            assert(i != line_out->a_id && i != line_out->b_id);
+            line_out->edge_type = ((m_use_quaternion ?
+                                    (m_quaternion * this->v_scaled_shared[i]).z()
+                                    : this->v_scaled_shared[i].z()) < slice_z) ? IntersectionLine::FacetEdgeType::Top : IntersectionLine::FacetEdgeType::Bottom;
+        }
+#endif
+        return FacetSliceType::Slicing;
+    }
+    return FacetSliceType::NoSlice;
+}
+
+static void slice_facet_at_zs(
+    const stl_triangle_vertex_indices                &indices,
+    const std::vector<Vec3f>                         &v_scaled_shared,
+    const int                                        *facet_neighbors,
+    const Eigen::Quaternion<float, Eigen::DontAlign> *quaternion,
+    std::vector<IntersectionLines>                   *lines,
+    boost::mutex                                     *lines_mutex,
+    const std::vector<float>                         &scaled_zs)
+{
+    stl_vertex                         vertices[3] { v_scaled_shared[indices(0)], v_scaled_shared[indices(1)], v_scaled_shared[indices(2)] };
+    if (quaternion)
+        for (int i = 0; i < 3; ++ i)
+            vertices[i] = *quaternion * vertices[i];
+
+    // find facet extents
+    const float min_z = fminf(vertices[0].z(), fminf(vertices[1].z(), vertices[2].z()));
+    const float max_z = fmaxf(vertices[0].z(), fmaxf(vertices[1].z(), vertices[2].z()));
+    
+    // find layer extents
+    auto min_layer = std::lower_bound(scaled_zs.begin(), scaled_zs.end(), min_z); // first layer whose slice_z is >= min_z
+    auto max_layer = std::upper_bound(min_layer, scaled_zs.end(), max_z); // first layer whose slice_z is > max_z
+    
+    for (auto it = min_layer; it != max_layer; ++ it) {
+        IntersectionLine il;
+        int idx_vertex_lowest = (vertices[1].z() == min_z) ? 1 : ((vertices[2].z() == min_z) ? 2 : 0);
+        if (slice_facet(*it, vertices, indices, facet_neighbors, idx_vertex_lowest, min_z == max_z, &il) == FacetSliceType::Slicing &&
+            il.edge_type != IntersectionLine::FacetEdgeType::Horizontal) {
+            // Ignore horizontal triangles. Any valid horizontal triangle must have a vertical triangle connected, otherwise the part has zero volume.
+            boost::lock_guard<boost::mutex> l(*lines_mutex);
+            (*lines)[it - scaled_zs.begin()].emplace_back(il);
+        }
+    }
+}
+
+#if 0
+//FIXME Should this go away? For valid meshes the function slice_facet() returns Slicing
+// and sets edges of vertical triangles to produce only a single edge per pair of neighbor faces.
+// So the following code makes only sense now to handle degenerate meshes with more than two faces
+// sharing a single edge.
+static inline void remove_tangent_edges(std::vector<IntersectionLine> &lines)
+{
+    std::vector<IntersectionLine*> by_vertex_pair;
+    by_vertex_pair.reserve(lines.size());
+    for (IntersectionLine& line : lines)
+        if (line.edge_type != IntersectionLine::FacetEdgeType::General && line.a_id != -1)
+            // This is a face edge. Check whether there is its neighbor stored in lines.
+            by_vertex_pair.emplace_back(&line);
+    auto edges_lower_sorted = [](const IntersectionLine *l1, const IntersectionLine *l2) {
+        // Sort vertices of l1, l2 lexicographically
+        int l1a = l1->a_id;
+        int l1b = l1->b_id;
+        int l2a = l2->a_id;
+        int l2b = l2->b_id;
+        if (l1a > l1b)
+            std::swap(l1a, l1b);
+        if (l2a > l2b)
+            std::swap(l2a, l2b);
+        // Lexicographical "lower" operator on lexicographically sorted vertices should bring equal edges together when sored.
+        return l1a < l2a || (l1a == l2a && l1b < l2b);
+    };
+    std::sort(by_vertex_pair.begin(), by_vertex_pair.end(), edges_lower_sorted);
+    for (auto line = by_vertex_pair.begin(); line != by_vertex_pair.end(); ++ line) {
+        IntersectionLine &l1 = **line;
+        if (! l1.skip()) {
+            // Iterate as long as line and line2 edges share the same end points.
+            for (auto line2 = line + 1; line2 != by_vertex_pair.end() && ! edges_lower_sorted(*line, *line2); ++ line2) {
+                // Lines must share the end points.
+                assert(! edges_lower_sorted(*line, *line2));
+                assert(! edges_lower_sorted(*line2, *line));
+                IntersectionLine &l2 = **line2;
+                if (l2.skip())
+                    continue;
+                if (l1.a_id == l2.a_id) {
+                    assert(l1.b_id == l2.b_id);
+                    l2.set_skip();
+                    // If they are both oriented upwards or downwards (like a 'V'),
+                    // then we can remove both edges from this layer since it won't 
+                    // affect the sliced shape.
+                    // If one of them is oriented upwards and the other is oriented
+                    // downwards, let's only keep one of them (it doesn't matter which
+                    // one since all 'top' lines were reversed at slicing).
+                    if (l1.edge_type == l2.edge_type) {
+                        l1.set_skip();
+                        break;
+                    }
+                } else {
+                    assert(l1.a_id == l2.b_id && l1.b_id == l2.a_id);
+                    // If this edge joins two horizontal facets, remove both of them.
+                    if (l1.edge_type == IntersectionLine::FacetEdgeType::Horizontal && l2.edge_type == IntersectionLine::FacetEdgeType::Horizontal) {
+                        l1.set_skip();
+                        l2.set_skip();
+                        break;
+                    }
+                }
+            }
+        }
+    }
+}
+#endif
+
+struct OpenPolyline {
+    OpenPolyline() = default;
+    OpenPolyline(const IntersectionReference &start, const IntersectionReference &end, Points &&points) : 
+        start(start), end(end), points(std::move(points)), consumed(false) { this->length = Slic3r::length(this->points); }
+    void reverse() {
+        std::swap(start, end);
+        std::reverse(points.begin(), points.end());
+    }
+    IntersectionReference   start;
+    IntersectionReference   end;
+    Points                  points;
+    double                  length;
+    bool                    consumed;
+};
+
+// called by TriangleMeshSlicer::make_loops() to connect sliced triangles into closed loops and open polylines by the triangle connectivity.
+// Only connects segments crossing triangles of the same orientation.
+static void chain_lines_by_triangle_connectivity(std::vector<IntersectionLine> &lines, Polygons &loops, std::vector<OpenPolyline> &open_polylines)
+{
+    // Build a map of lines by edge_a_id and a_id.
+    std::vector<IntersectionLine*> by_edge_a_id;
+    std::vector<IntersectionLine*> by_a_id;
+    by_edge_a_id.reserve(lines.size());
+    by_a_id.reserve(lines.size());
+    for (IntersectionLine &line : lines) {
+        if (! line.skip()) {
+            if (line.edge_a_id != -1)
+                by_edge_a_id.emplace_back(&line);
+            if (line.a_id != -1)
+                by_a_id.emplace_back(&line);
+        }
+    }
+    auto by_edge_lower = [](const IntersectionLine* il1, const IntersectionLine *il2) { return il1->edge_a_id < il2->edge_a_id; };
+    auto by_vertex_lower = [](const IntersectionLine* il1, const IntersectionLine *il2) { return il1->a_id < il2->a_id; };
+    std::sort(by_edge_a_id.begin(), by_edge_a_id.end(), by_edge_lower);
+    std::sort(by_a_id.begin(), by_a_id.end(), by_vertex_lower);
+    // Chain the segments with a greedy algorithm, collect the loops and unclosed polylines.
+    IntersectionLines::iterator it_line_seed = lines.begin();
+    for (;;) {
+        // take first spare line and start a new loop
+        IntersectionLine *first_line = nullptr;
+        for (; it_line_seed != lines.end(); ++ it_line_seed)
+            if (it_line_seed->is_seed_candidate()) {
+            //if (! it_line_seed->skip()) {
+                first_line = &(*it_line_seed ++);
+                break;
+            }
+        if (first_line == nullptr)
+            break;
+        first_line->set_skip();
+        Points loop_pts;
+        loop_pts.emplace_back(first_line->a);
+        IntersectionLine *last_line = first_line;
+        
+        /*
+        printf("first_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", 
+            first_line->edge_a_id, first_line->edge_b_id, first_line->a_id, first_line->b_id,
+            first_line->a.x, first_line->a.y, first_line->b.x, first_line->b.y);
+        */
+        
+        IntersectionLine key;
+        for (;;) {
+            // find a line starting where last one finishes
+            IntersectionLine* next_line = nullptr;
+            if (last_line->edge_b_id != -1) {
+                key.edge_a_id = last_line->edge_b_id;
+                auto it_begin = std::lower_bound(by_edge_a_id.begin(), by_edge_a_id.end(), &key, by_edge_lower);
+                if (it_begin != by_edge_a_id.end()) {
+                    auto it_end = std::upper_bound(it_begin, by_edge_a_id.end(), &key, by_edge_lower);
+                    for (auto it_line = it_begin; it_line != it_end; ++ it_line)
+                        if (! (*it_line)->skip()) {
+                            next_line = *it_line;
+                            break;
+                        }
+                }
+            }
+            if (next_line == nullptr && last_line->b_id != -1) {
+                key.a_id = last_line->b_id;
+                auto it_begin = std::lower_bound(by_a_id.begin(), by_a_id.end(), &key, by_vertex_lower);
+                if (it_begin != by_a_id.end()) {
+                    auto it_end = std::upper_bound(it_begin, by_a_id.end(), &key, by_vertex_lower);
+                    for (auto it_line = it_begin; it_line != it_end; ++ it_line)
+                        if (! (*it_line)->skip()) {
+                            next_line = *it_line;
+                            break;
+                        }
+                }
+            }
+            if (next_line == nullptr) {
+                // Check whether we closed this loop.
+                if ((first_line->edge_a_id != -1 && first_line->edge_a_id == last_line->edge_b_id) || 
+                    (first_line->a_id      != -1 && first_line->a_id      == last_line->b_id)) {
+                    // The current loop is complete. Add it to the output.
+                    loops.emplace_back(std::move(loop_pts));
+                    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
+                    printf("  Discovered %s polygon of %d points\n", (p.is_counter_clockwise() ? "ccw" : "cw"), (int)p.points.size());
+                    #endif
+                } else {
+                    // This is an open polyline. Add it to the list of open polylines. These open polylines will processed later.
+                    loop_pts.emplace_back(last_line->b);
+                    open_polylines.emplace_back(OpenPolyline(
+                        IntersectionReference(first_line->a_id, first_line->edge_a_id), 
+                        IntersectionReference(last_line->b_id, last_line->edge_b_id), std::move(loop_pts)));
+                }
+                break;
+            }
+            /*
+            printf("next_line edge_a_id = %d, edge_b_id = %d, a_id = %d, b_id = %d, a = %d,%d, b = %d,%d\n", 
+                next_line->edge_a_id, next_line->edge_b_id, next_line->a_id, next_line->b_id,
+                next_line->a.x, next_line->a.y, next_line->b.x, next_line->b.y);
+            */
+            loop_pts.emplace_back(next_line->a);
+            last_line = next_line;
+            next_line->set_skip();
+        }
+    }
+}
+
+std::vector<OpenPolyline*> open_polylines_sorted(std::vector<OpenPolyline> &open_polylines, bool update_lengths)
+{
+    std::vector<OpenPolyline*> out;
+    out.reserve(open_polylines.size());
+    for (OpenPolyline &opl : open_polylines)
+        if (! opl.consumed) {
+            if (update_lengths)
+                opl.length = Slic3r::length(opl.points);
+            out.emplace_back(&opl);
+        }
+    std::sort(out.begin(), out.end(), [](const OpenPolyline *lhs, const OpenPolyline *rhs){ return lhs->length > rhs->length; });
+    return out;
+}
+
+// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices.
+// Depending on "try_connect_reversed", it may or may not connect segments crossing triangles of opposite orientation.
+static void chain_open_polylines_exact(std::vector<OpenPolyline> &open_polylines, Polygons &loops, bool try_connect_reversed)
+{
+    // Store the end points of open_polylines into vectors sorted
+    struct OpenPolylineEnd {
+        OpenPolylineEnd(OpenPolyline *polyline, bool start) : polyline(polyline), start(start) {}
+        OpenPolyline    *polyline;
+        // Is it the start or end point?
+        bool             start;
+        const IntersectionReference& ipref() const { return start ? polyline->start : polyline->end; }
+        // Return a unique ID for the intersection point.
+        // Return a positive id for a point, or a negative id for an edge.
+        int id() const { const IntersectionReference &r = ipref(); return (r.point_id >= 0) ? r.point_id : - r.edge_id; }
+        bool operator==(const OpenPolylineEnd &rhs) const { return this->polyline == rhs.polyline && this->start == rhs.start; }
+    };
+    auto by_id_lower = [](const OpenPolylineEnd &ope1, const OpenPolylineEnd &ope2) { return ope1.id() < ope2.id(); };
+    std::vector<OpenPolylineEnd> by_id;
+    by_id.reserve(2 * open_polylines.size());
+    for (OpenPolyline &opl : open_polylines) {
+        if (opl.start.point_id != -1 || opl.start.edge_id != -1)
+            by_id.emplace_back(OpenPolylineEnd(&opl, true));
+        if (try_connect_reversed && (opl.end.point_id != -1 || opl.end.edge_id != -1))
+            by_id.emplace_back(OpenPolylineEnd(&opl, false));
+    }
+    std::sort(by_id.begin(), by_id.end(), by_id_lower);
+    // Find an iterator to by_id_lower for the particular end of OpenPolyline (by comparing the OpenPolyline pointer and the start attribute).
+    auto find_polyline_end = [&by_id, by_id_lower](const OpenPolylineEnd &end) -> std::vector<OpenPolylineEnd>::iterator {
+        for (auto it = std::lower_bound(by_id.begin(), by_id.end(), end, by_id_lower);
+                  it != by_id.end() && it->id() == end.id(); ++ it)
+            if (*it == end)
+                return it;
+        return by_id.end();
+    };
+    // Try to connect the loops.
+    std::vector<OpenPolyline*> sorted_by_length = open_polylines_sorted(open_polylines, false);
+    for (OpenPolyline *opl : sorted_by_length) {
+        if (opl->consumed)
+            continue;
+        opl->consumed = true;
+        OpenPolylineEnd end(opl, false);
+        for (;;) {
+            // find a line starting where last one finishes
+            auto it_next_start = std::lower_bound(by_id.begin(), by_id.end(), end, by_id_lower);
+            for (; it_next_start != by_id.end() && it_next_start->id() == end.id(); ++ it_next_start)
+                if (! it_next_start->polyline->consumed)
+                    goto found;
+            // The current loop could not be closed. Unmark the segment.
+            opl->consumed = false;
+            break;
+        found:
+            // Attach this polyline to the end of the initial polyline.
+            if (it_next_start->start) {
+                auto it = it_next_start->polyline->points.begin();
+                std::copy(++ it, it_next_start->polyline->points.end(), back_inserter(opl->points));
+            } else {
+                auto it = it_next_start->polyline->points.rbegin();
+                std::copy(++ it, it_next_start->polyline->points.rend(), back_inserter(opl->points));
+            }
+            opl->length += it_next_start->polyline->length;
+            // Mark the next polyline as consumed.
+            it_next_start->polyline->points.clear();
+            it_next_start->polyline->length = 0.;
+            it_next_start->polyline->consumed = true;
+            if (try_connect_reversed) {
+                // Running in a mode, where the polylines may be connected by mixing their orientations.
+                // Update the end point lookup structure after the end point of the current polyline was extended.
+                auto it_end      = find_polyline_end(end);
+                auto it_next_end = find_polyline_end(OpenPolylineEnd(it_next_start->polyline, !it_next_start->start));
+                // Swap the end points of the current and next polyline, but keep the polyline ptr and the start flag.
+                std::swap(opl->end, it_next_end->start ? it_next_end->polyline->start : it_next_end->polyline->end);
+                // Swap the positions of OpenPolylineEnd structures in the sorted array to match their respective end point positions.
+                std::swap(*it_end, *it_next_end);
+            }
+            // Check whether we closed this loop.
+            if ((opl->start.edge_id  != -1 && opl->start.edge_id  == opl->end.edge_id) ||
+                (opl->start.point_id != -1 && opl->start.point_id == opl->end.point_id)) {
+                // The current loop is complete. Add it to the output.
+                //assert(opl->points.front().point_id == opl->points.back().point_id);
+                //assert(opl->points.front().edge_id  == opl->points.back().edge_id);
+                // Remove the duplicate last point.
+                opl->points.pop_back();
+                if (opl->points.size() >= 3) {
+                    if (try_connect_reversed && area(opl->points) < 0)
+                        // The closed polygon is patched from pieces with messed up orientation, therefore
+                        // the orientation of the patched up polygon is not known.
+                        // Orient the patched up polygons CCW. This heuristic may close some holes and cavities.
+                        std::reverse(opl->points.begin(), opl->points.end());
+                    loops.emplace_back(std::move(opl->points));
+                }
+                opl->points.clear();
+                break;
+            }
+            // Continue with the current loop.
+        }
+    }
+}
+
+// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices, 
+// possibly closing small gaps.
+// Depending on "try_connect_reversed", it may or may not connect segments crossing triangles of opposite orientation.
+static void chain_open_polylines_close_gaps(std::vector<OpenPolyline> &open_polylines, Polygons &loops, double max_gap, bool try_connect_reversed)
+{
+    const coord_t max_gap_scaled = (coord_t)scale_(max_gap);
+
+    // Sort the open polylines by their length, so the new loops will be seeded from longer chains.
+    // Update the polyline lengths, return only not yet consumed polylines.
+    std::vector<OpenPolyline*> sorted_by_length = open_polylines_sorted(open_polylines, true);
+
+    // Store the end points of open_polylines into ClosestPointInRadiusLookup<OpenPolylineEnd>.
+    struct OpenPolylineEnd {
+        OpenPolylineEnd(OpenPolyline *polyline, bool start) : polyline(polyline), start(start) {}
+        OpenPolyline    *polyline;
+        // Is it the start or end point?
+        bool             start;
+        const Point&     point() const { return start ? polyline->points.front() : polyline->points.back(); }
+        bool operator==(const OpenPolylineEnd &rhs) const { return this->polyline == rhs.polyline && this->start == rhs.start; }
+    };
+    struct OpenPolylineEndAccessor {
+        const Point* operator()(const OpenPolylineEnd &pt) const { return pt.polyline->consumed ? nullptr : &pt.point(); }
+    };
+    typedef ClosestPointInRadiusLookup<OpenPolylineEnd, OpenPolylineEndAccessor> ClosestPointLookupType;
+    ClosestPointLookupType closest_end_point_lookup(max_gap_scaled);
+    for (OpenPolyline *opl : sorted_by_length) {
+        closest_end_point_lookup.insert(OpenPolylineEnd(opl, true));
+        if (try_connect_reversed)
+            closest_end_point_lookup.insert(OpenPolylineEnd(opl, false));
+    }
+    // Try to connect the loops.
+    for (OpenPolyline *opl : sorted_by_length) {
+        if (opl->consumed)
+            continue;
+        OpenPolylineEnd end(opl, false);
+        if (try_connect_reversed)
+            // The end point of this polyline will be modified, thus the following entry will become invalid. Remove it.
+            closest_end_point_lookup.erase(end);
+        opl->consumed = true;
+        size_t n_segments_joined = 1;
+        for (;;) {
+            // Find a line starting where last one finishes, only return non-consumed open polylines (OpenPolylineEndAccessor returns null for consumed).
+            std::pair<const OpenPolylineEnd*, double> next_start_and_dist = closest_end_point_lookup.find(end.point());
+            const OpenPolylineEnd *next_start = next_start_and_dist.first;
+            // Check whether we closed this loop.
+            double current_loop_closing_distance2 = (opl->points.back() - opl->points.front()).cast<double>().squaredNorm();
+            bool   loop_closed = current_loop_closing_distance2 < coordf_t(max_gap_scaled) * coordf_t(max_gap_scaled);
+            if (next_start != nullptr && loop_closed && current_loop_closing_distance2 < next_start_and_dist.second) {
+                // Heuristics to decide, whether to close the loop, or connect another polyline.
+                // One should avoid closing loops shorter than max_gap_scaled.
+                loop_closed = sqrt(current_loop_closing_distance2) < 0.3 * length(opl->points);
+            }
+            if (loop_closed) {
+                // Remove the start point of the current polyline from the lookup.
+                // Mark the current segment as not consumed, otherwise the closest_end_point_lookup.erase() would fail.
+                opl->consumed = false;
+                closest_end_point_lookup.erase(OpenPolylineEnd(opl, true));
+                if (current_loop_closing_distance2 == 0.) {
+                    // Remove the duplicate last point.
+                    opl->points.pop_back();
+                } else {
+                    // The end points are different, keep both of them.
+                }
+                if (opl->points.size() >= 3) {
+                    if (try_connect_reversed && n_segments_joined > 1 && area(opl->points) < 0)
+                        // The closed polygon is patched from pieces with messed up orientation, therefore
+                        // the orientation of the patched up polygon is not known.
+                        // Orient the patched up polygons CCW. This heuristic may close some holes and cavities.
+                        std::reverse(opl->points.begin(), opl->points.end());
+                    loops.emplace_back(std::move(opl->points));
+                }
+                opl->points.clear();
+                opl->consumed = true;
+                break;
+            }
+            if (next_start == nullptr) {
+                // The current loop could not be closed. Unmark the segment.
+                opl->consumed = false;
+                if (try_connect_reversed)
+                    // Re-insert the end point.
+                    closest_end_point_lookup.insert(OpenPolylineEnd(opl, false));
+                break;
+            }
+            // Attach this polyline to the end of the initial polyline.
+            if (next_start->start) {
+                auto it = next_start->polyline->points.begin();
+                if (*it == opl->points.back())
+                    ++ it;
+                std::copy(it, next_start->polyline->points.end(), back_inserter(opl->points));
+            } else {
+                auto it = next_start->polyline->points.rbegin();
+                if (*it == opl->points.back())
+                    ++ it;
+                std::copy(it, next_start->polyline->points.rend(), back_inserter(opl->points));
+            }
+            ++ n_segments_joined;
+            // Remove the end points of the consumed polyline segment from the lookup.
+            OpenPolyline *opl2 = next_start->polyline;
+            closest_end_point_lookup.erase(OpenPolylineEnd(opl2, true));
+            if (try_connect_reversed)
+                closest_end_point_lookup.erase(OpenPolylineEnd(opl2, false));
+            opl2->points.clear();
+            opl2->consumed = true;
+            // Continue with the current loop.
+        }
+    }
+}
+
+static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
+{
+#if 0
+//FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
+//#ifdef _DEBUG
+    for (const Line &l : lines)
+        assert(l.a != l.b);
+#endif /* _DEBUG */
+
+    // There should be no tangent edges, as the horizontal triangles are ignored and if two triangles touch at a cutting plane,
+    // only the bottom triangle is considered to be cutting the plane.
+//    remove_tangent_edges(lines);
+
+#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
+        BoundingBox bbox_svg;
+        {
+            static int iRun = 0;
+            for (const Line &line : lines) {
+                bbox_svg.merge(line.a);
+                bbox_svg.merge(line.b);
+            }
+            SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-raw_lines-%d.svg", iRun ++).c_str(), bbox_svg);
+            for (const Line &line : lines)
+                svg.draw(line);
+            svg.Close();
+        }
+#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+
+    std::vector<OpenPolyline> open_polylines;
+    chain_lines_by_triangle_connectivity(lines, *loops, open_polylines);
+
+#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
+        {
+            static int iRun = 0;
+            SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines-%d.svg", iRun ++).c_str(), bbox_svg);
+            svg.draw(union_ex(*loops));
+            for (const OpenPolyline &pl : open_polylines)
+                svg.draw(Polyline(pl.points), "red");
+            svg.Close();
+        }
+#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+
+    // Now process the open polylines.
+    // Do it in two rounds, first try to connect in the same direction only,
+    // then try to connect the open polylines in reversed order as well.
+    chain_open_polylines_exact(open_polylines, *loops, false);
+    chain_open_polylines_exact(open_polylines, *loops, true);
+
+#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
+    {
+        static int iRun = 0;
+        SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines2-%d.svg", iRun++).c_str(), bbox_svg);
+        svg.draw(union_ex(*loops));
+        for (const OpenPolyline &pl : open_polylines) {
+            if (pl.points.empty())
+                continue;
+            svg.draw(Polyline(pl.points), "red");
+            svg.draw(pl.points.front(), "blue");
+            svg.draw(pl.points.back(), "blue");
+        }
+        svg.Close();
+    }
+#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+
+    // Try to close gaps.
+    // Do it in two rounds, first try to connect in the same direction only,
+    // then try to connect the open polylines in reversed order as well.
+#if 0
+    for (double max_gap : { EPSILON, 0.001, 0.1, 1., 2. }) {
+        chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, false);
+        chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, true);
+    }
+#else
+    const double max_gap = 2.; //mm
+    chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, false);
+    chain_open_polylines_close_gaps(open_polylines, *loops, max_gap, true);
+#endif
+
+#ifdef SLIC3R_DEBUG_SLICE_PROCESSING
+    {
+        static int iRun = 0;
+        SVG svg(debug_out_path("TriangleMeshSlicer_make_loops-polylines-final-%d.svg", iRun++).c_str(), bbox_svg);
+        svg.draw(union_ex(*loops));
+        for (const OpenPolyline &pl : open_polylines) {
+            if (pl.points.empty())
+                continue;
+            svg.draw(Polyline(pl.points), "red");
+            svg.draw(pl.points.front(), "blue");
+            svg.draw(pl.points.back(), "blue");
+        }
+        svg.Close();
+    }
+#endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+}
+
+// Used to cut the mesh into two halves.
+static ExPolygons make_expolygons_simple(std::vector<IntersectionLine> &lines)
+{
+    ExPolygons slices;
+    Polygons holes;
+
+    {
+        Polygons loops;
+        make_loops(lines, &loops);        
+        for (Polygon &loop : loops)
+            if (loop.area() >= 0.)
+                slices.emplace_back(std::move(loop));
+            else
+                holes.emplace_back(std::move(loop));
+    }
+
+    // If there are holes, then there should also be outer contours.
+    assert(holes.empty() || ! slices.empty());
+    if (! slices.empty())
+    {
+        // Assign holes to outer contours.
+        for (Polygon &hole : holes) {
+            // Find an outer contour to a hole.
+            int     slice_idx            = -1;
+            double  current_contour_area = std::numeric_limits<double>::max();
+            for (ExPolygon &slice : slices)
+                if (slice.contour.contains(hole.points.front())) {
+                    double area = slice.contour.area();
+                    if (area < current_contour_area) {
+                        slice_idx = &slice - slices.data();
+                        current_contour_area = area;
+                    }
+                }
+            // assert(slice_idx != -1);
+            if (slice_idx == -1)
+                // Ignore this hole.
+                continue;
+            assert(current_contour_area < std::numeric_limits<double>::max() && current_contour_area >= -hole.area());
+            slices[slice_idx].holes.emplace_back(std::move(hole));
+        }
+
+#if 0
+        // If the input mesh is not valid, the holes may intersect with the external contour.
+        // Rather subtract them from the outer contour.
+        Polygons poly;
+        for (auto it_slice = slices->begin(); it_slice != slices->end(); ++ it_slice) {
+            if (it_slice->holes.empty()) {
+                poly.emplace_back(std::move(it_slice->contour));
+            } else {
+                Polygons contours;
+                contours.emplace_back(std::move(it_slice->contour));
+                for (auto it = it_slice->holes.begin(); it != it_slice->holes.end(); ++ it)
+                    it->reverse();
+                polygons_append(poly, diff(contours, it_slice->holes));
+            }
+        }
+        // If the input mesh is not valid, the input contours may intersect.
+        *slices = union_ex(poly);
+#endif
+
+#if 0
+        // If the input mesh is not valid, the holes may intersect with the external contour.
+        // Rather subtract them from the outer contour.
+        ExPolygons poly;
+        for (auto it_slice = slices->begin(); it_slice != slices->end(); ++ it_slice) {
+            Polygons contours;
+            contours.emplace_back(std::move(it_slice->contour));
+            for (auto it = it_slice->holes.begin(); it != it_slice->holes.end(); ++ it)
+                it->reverse();
+            expolygons_append(poly, diff_ex(contours, it_slice->holes));
+        }
+        // If the input mesh is not valid, the input contours may intersect.
+        *slices = std::move(poly);
+#endif
+    }
+
+    return slices;
+}
+
+static void make_expolygons(const Polygons &loops, const float closing_radius, const float extra_offset, ExPolygons* slices)
+{
+    /*
+        Input loops are not suitable for evenodd nor nonzero fill types, as we might get
+        two consecutive concentric loops having the same winding order - and we have to 
+        respect such order. In that case, evenodd would create wrong inversions, and nonzero
+        would ignore holes inside two concentric contours.
+        So we're ordering loops and collapse consecutive concentric loops having the same 
+        winding order.
+        TODO: find a faster algorithm for this, maybe with some sort of binary search.
+        If we computed a "nesting tree" we could also just remove the consecutive loops
+        having the same winding order, and remove the extra one(s) so that we could just
+        supply everything to offset() instead of performing several union/diff calls.
+    
+        we sort by area assuming that the outermost loops have larger area;
+        the previous sorting method, based on $b->contains($a->[0]), failed to nest
+        loops correctly in some edge cases when original model had overlapping facets
+    */
+
+    /* The following lines are commented out because they can generate wrong polygons,
+       see for example issue #661 */
+
+    //std::vector<double> area;
+    //std::vector<size_t> sorted_area;  // vector of indices
+    //for (Polygons::const_iterator loop = loops.begin(); loop != loops.end(); ++ loop) {
+    //    area.emplace_back(loop->area());
+    //    sorted_area.emplace_back(loop - loops.begin());
+    //}
+    //
+    //// outer first
+    //std::sort(sorted_area.begin(), sorted_area.end(),
+    //    [&area](size_t a, size_t b) { return std::abs(area[a]) > std::abs(area[b]); });
+
+    //// we don't perform a safety offset now because it might reverse cw loops
+    //Polygons p_slices;
+    //for (std::vector<size_t>::const_iterator loop_idx = sorted_area.begin(); loop_idx != sorted_area.end(); ++ loop_idx) {
+    //    /* we rely on the already computed area to determine the winding order
+    //       of the loops, since the Orientation() function provided by Clipper
+    //       would do the same, thus repeating the calculation */
+    //    Polygons::const_iterator loop = loops.begin() + *loop_idx;
+    //    if (area[*loop_idx] > +EPSILON)
+    //        p_slices.emplace_back(*loop);
+    //    else if (area[*loop_idx] < -EPSILON)
+    //        //FIXME This is arbitrary and possibly very slow.
+    //        // If the hole is inside a polygon, then there is no need to diff.
+    //        // If the hole intersects a polygon boundary, then diff it, but then
+    //        // there is no guarantee of an ordering of the loops.
+    //        // Maybe we can test for the intersection before running the expensive diff algorithm?
+    //        p_slices = diff(p_slices, *loop);
+    //}
+
+    // Perform a safety offset to merge very close facets (TODO: find test case for this)
+    // 0.0499 comes from https://github.com/slic3r/Slic3r/issues/959
+//    double safety_offset = scale_(0.0499);
+    // 0.0001 is set to satisfy GH #520, #1029, #1364
+    assert(closing_radius >= 0);
+    assert(extra_offset >= 0);
+    double offset_out = + scale_(closing_radius + extra_offset);
+    double offset_in  = - scale_(closing_radius);
+
+    /* The following line is commented out because it can generate wrong polygons,
+       see for example issue #661 */
+    //ExPolygons ex_slices = offset2_ex(p_slices, +safety_offset, -safety_offset);
+    
+    #ifdef SLIC3R_TRIANGLEMESH_DEBUG
+    size_t holes_count = 0;
+    for (ExPolygons::const_iterator e = ex_slices.begin(); e != ex_slices.end(); ++ e)
+        holes_count += e->holes.size();
+    printf("%zu surface(s) having %zu holes detected from %zu polylines\n",
+        ex_slices.size(), holes_count, loops.size());
+    #endif
+    
+    // append to the supplied collection
+    expolygons_append(*slices,
+        offset_out > 0 && offset_in < 0 ? offset2_ex(union_ex(loops), offset_out, offset_in) :
+        offset_out > 0 ? offset_ex(union_ex(loops), offset_out) :
+        offset_in  < 0 ? offset_ex(union_ex(loops), offset_in) :
+        union_ex(loops));
+}
+
+/*
+static void make_expolygons(std::vector<IntersectionLine> &lines, const float closing_radius, ExPolygons* slices)
+{
+    Polygons pp;
+    make_loops(lines, &pp);
+    Slic3r::make_expolygons(pp, closing_radius, 0.f, slices);
+}
+*/
+
+void TriangleMeshSlicer::init(const TriangleMesh *mesh, throw_on_cancel_callback_type throw_on_cancel)
+{
+    if (! mesh->has_shared_vertices())
+        throw Slic3r::InvalidArgument("TriangleMeshSlicer was passed a mesh without shared vertices.");
+
+    this->init(&mesh->its, throw_on_cancel);
+}
+
+void TriangleMeshSlicer::init(const indexed_triangle_set *its, throw_on_cancel_callback_type throw_on_cancel)
+{
+    m_its = its;
+    facets_edges = create_face_neighbors_index(*its, throw_on_cancel);
+    v_scaled_shared.assign(its->vertices.size(), stl_vertex());
+    for (size_t i = 0; i < v_scaled_shared.size(); ++ i)
+        this->v_scaled_shared[i] = its->vertices[i] / float(SCALING_FACTOR);
+}
+
+void TriangleMeshSlicer::set_up_direction(const Vec3f& up)
+{
+    m_quaternion.setFromTwoVectors(up, Vec3f::UnitZ());
+    m_use_quaternion = true;
+}
+
+void TriangleMeshSlicer::slice(
+    const std::vector<float>         &z,
+    const MeshSlicingParams          &params,
+    std::vector<Polygons>           *layers,
+    throw_on_cancel_callback_type    throw_on_cancel) const
+{
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice";
+
+    /*
+       This method gets called with a list of unscaled Z coordinates and outputs
+       a vector pointer having the same number of items as the original list.
+       Each item is a vector of polygons created by slicing our mesh at the 
+       given heights.
+       
+       This method should basically combine the behavior of the existing
+       Perl methods defined in lib/Slic3r/TriangleMesh.pm:
+       
+       - analyze(): this creates the 'facets_edges' and the 'edges_facets'
+            tables (we don't need the 'edges' table)
+       
+       - slice_facet(): this has to be done for each facet. It generates 
+            intersection lines with each plane identified by the Z list.
+            The get_layer_range() binary search used to identify the Z range
+            of the facet is already ported to C++ (see Object.xsp)
+       
+       - make_loops(): this has to be done for each layer. It creates polygons
+            from the lines generated by the previous step.
+        
+        At the end, we free the tables generated by analyze() as we don't 
+        need them anymore.
+        
+        NOTE: this method accepts a vector of floats because the mesh coordinate
+        type is float.
+    */
+    
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice_facet_at_zs";
+    std::vector<IntersectionLines> lines(z.size());
+    {
+        std::vector<float> scaled_z(z);
+        for (float &z : scaled_z)
+            z = scaled<float>(z);
+        boost::mutex lines_mutex;
+        tbb::parallel_for(
+            tbb::blocked_range<int>(0, int(m_its->indices.size())),
+            [&lines, &lines_mutex, &scaled_z, throw_on_cancel, this](const tbb::blocked_range<int>& range) {
+                const Eigen::Quaternion<float, Eigen::DontAlign> *rotation = m_use_quaternion ? &m_quaternion : nullptr;
+                for (int facet_idx = range.begin(); facet_idx < range.end(); ++ facet_idx) {
+                    if ((facet_idx & 0x0ffff) == 0)
+                        throw_on_cancel();
+                    slice_facet_at_zs(m_its->indices[facet_idx], this->v_scaled_shared, this->facets_edges.data() + facet_idx * 3, rotation, &lines, &lines_mutex, scaled_z);
+                }
+            }
+        );
+    }
+    throw_on_cancel();
+
+    // v_scaled_shared could be freed here
+    
+    // build loops
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::_make_loops_do";
+    layers->resize(z.size());
+    tbb::parallel_for(
+        tbb::blocked_range<size_t>(0, z.size()),
+        [&lines, &layers, &params, throw_on_cancel](const tbb::blocked_range<size_t>& range) {
+            for (size_t line_idx = range.begin(); line_idx < range.end(); ++ line_idx) {
+                if ((line_idx & 0x0ffff) == 0)
+                    throw_on_cancel();
+
+                Polygons   &polygons  = (*layers)[line_idx];
+                make_loops(lines[line_idx], &polygons);
+
+                auto this_mode = line_idx < params.slicing_mode_normal_below_layer ? params.mode_below : params.mode;
+                if (! polygons.empty()) {
+                    if (this_mode == SlicingMode::Positive) {
+                        // Reorient all loops to be CCW.
+                        for (Polygon& p : polygons)
+                            p.make_counter_clockwise();
+                    } else if (this_mode == SlicingMode::PositiveLargestContour) {
+                        // Keep just the largest polygon, make it CCW.
+                        double   max_area = 0.;
+                        Polygon* max_area_polygon = nullptr;
+                        for (Polygon& p : polygons) {
+                            double a = p.area();
+                            if (std::abs(a) > std::abs(max_area)) {
+                                max_area = a;
+                                max_area_polygon = &p;
+                            }
+                        }
+                        assert(max_area_polygon != nullptr);
+                        if (max_area < 0.)
+                            max_area_polygon->reverse();
+                        Polygon p(std::move(*max_area_polygon));
+                        polygons.clear();
+                        polygons.emplace_back(std::move(p));
+                    }
+                }
+            }
+        }
+    );
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice finished";
+
+#ifdef SLIC3R_DEBUG
+    {
+        static int iRun = 0;
+        for (size_t i = 0; i < z.size(); ++ i) {
+            Polygons  &polygons   = (*layers)[i];
+            ExPolygons expolygons = union_ex(polygons, true);
+            SVG::export_expolygons(debug_out_path("slice_%d_%d.svg", iRun, i).c_str(), expolygons);
+            {
+                BoundingBox bbox;
+                for (const IntersectionLine &l : lines[i]) {
+                    bbox.merge(l.a);
+                    bbox.merge(l.b);
+                }
+                SVG svg(debug_out_path("slice_loops_%d_%d.svg", iRun, i).c_str(), bbox);
+                svg.draw(expolygons);
+                for (const IntersectionLine &l : lines[i])
+                    svg.draw(l, "red", 0);
+                svg.draw_outline(expolygons, "black", "blue", 0);
+                svg.Close();
+            }
+#if 0
+//FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
+            for (Polygon &poly : polygons) {
+                for (size_t i = 1; i < poly.points.size(); ++ i)
+                    assert(poly.points[i-1] != poly.points[i]);
+                assert(poly.points.front() != poly.points.back());
+            }
+#endif
+        }
+        ++ iRun;
+    }
+#endif
+}
+
+void TriangleMeshSlicer::slice(
+    // Where to slice.
+    const std::vector<float>        &z,
+    const MeshSlicingParamsExtended &params,
+    std::vector<ExPolygons>        *layers,
+    throw_on_cancel_callback_type   throw_on_cancel) const
+{
+    std::vector<Polygons> layers_p;
+    {
+        MeshSlicingParams slicing_params(params);
+        if (params.mode == SlicingMode::PositiveLargestContour)
+            slicing_params.mode = SlicingMode::Positive;
+        if (params.mode_below == SlicingMode::PositiveLargestContour)
+            slicing_params.mode_below = SlicingMode::Positive;
+        this->slice(z, slicing_params, &layers_p, throw_on_cancel);
+    }
+    
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - start";
+    layers->resize(z.size());
+    tbb::parallel_for(
+        tbb::blocked_range<size_t>(0, z.size()),
+        [&layers_p, &params, layers, throw_on_cancel]
+        (const tbb::blocked_range<size_t>& range) {
+            for (size_t layer_id = range.begin(); layer_id < range.end(); ++ layer_id) {
+#ifdef SLIC3R_TRIANGLEMESH_DEBUG
+                printf("Layer %zu (slice_z = %.2f):\n", layer_id, z[layer_id]);
+#endif
+                throw_on_cancel();
+                ExPolygons &expolygons = (*layers)[layer_id];
+                Slic3r::make_expolygons(layers_p[layer_id], params.closing_radius, params.extra_offset, &expolygons);
+                //FIXME simplify
+                const auto this_mode = layer_id < params.slicing_mode_normal_below_layer ? params.mode_below : params.mode;
+                if (this_mode == SlicingMode::PositiveLargestContour)
+                    keep_largest_contour_only(expolygons);
+            }
+        });
+    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - end";
+}
+
+static void triangulate_slice(indexed_triangle_set &its, IntersectionLines &lines, const std::vector<int> &slice_vertices, float z)
+{
+//    BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - triangulating the cut";
+
+    ExPolygons section   = make_expolygons_simple(lines);
+    Pointf3s   triangles = triangulate_expolygons_3d(section, z, true);
+    std::vector<std::pair<Vec2f, int>> map_vertex_to_index;
+    map_vertex_to_index.reserve(slice_vertices.size());
+    for (int i : slice_vertices)
+        map_vertex_to_index.emplace_back(to_2d(its.vertices[i]), i);
+    std::sort(map_vertex_to_index.begin(), map_vertex_to_index.end(), 
+        [](const std::pair<Vec2f, int> &l, const std::pair<Vec2f, int> &r) { 
+            return l.first.x() < r.first.x() || (l.first.x() == r.first.x() && l.first.y() < r.first.y()); });
+    size_t idx_vertex_new_first = its.vertices.size();
+    for (size_t i = 0; i < triangles.size(); ) {
+        stl_triangle_vertex_indices facet;
+        for (size_t j = 0; j < 3; ++ j) {
+            Vec3f v = triangles[i++].cast<float>();
+            auto it = lower_bound_by_predicate(map_vertex_to_index.begin(), map_vertex_to_index.end(), 
+                [&v](const std::pair<Vec2f, int> &l) { return l.first.x() < v.x() || (l.first.x() == v.x() && l.first.y() < v.y()); });
+            int   idx = -1;
+            if (it != map_vertex_to_index.end() && it->first.x() == v.x() && it->first.y() == v.y())
+                idx = it->second;
+            else {
+                // Try to find the vertex in the list of newly added vertices. Those vertices are not matched on the cut and they shall be rare.
+                for (size_t k = idx_vertex_new_first; k < its.vertices.size(); ++ k)
+                    if (its.vertices[k] == v) {
+                        idx = int(k);
+                        break;
+                    }
+                if (idx == -1) {
+                    idx = int(its.vertices.size());
+                    its.vertices.emplace_back(v);
+                }
+            }
+            facet(j) = idx;
+        }
+        its.indices.emplace_back(facet);
+    }
+
+    its_compactify_vertices(its);
+    its_remove_degenerate_faces(its);
+}
+
+void TriangleMeshSlicer::cut(float z, indexed_triangle_set *upper, indexed_triangle_set *lower) const
+{
+    assert(upper || lower);
+    if (upper == nullptr && lower == nullptr)
+        return;
+
+    if (upper) {
+        upper->clear();
+        upper->vertices = m_its->vertices;
+        upper->indices.reserve(m_its->indices.size());
+    }
+    if (lower) {
+        lower->clear();
+        lower->vertices = m_its->vertices;
+        lower->indices.reserve(m_its->indices.size());
+    }
+
+    BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - slicing object";
+    const auto scaled_z = scaled<float>(z);
+
+    // To triangulate the caps after slicing.
+    IntersectionLines upper_lines, lower_lines;
+    std::vector<int>  upper_slice_vertices, lower_slice_vertices;
+
+    for (int facet_idx = 0; facet_idx < int(m_its->indices.size()); ++ facet_idx) {
+        const stl_triangle_vertex_indices &facet = m_its->indices[facet_idx];
+        Vec3f vertices[3] { m_its->vertices[facet(0)], m_its->vertices[facet(1)], m_its->vertices[facet(2)] };
+        stl_vertex vertices_scaled[3]{ this->v_scaled_shared[facet[0]], this->v_scaled_shared[facet[1]], this->v_scaled_shared[facet[2]] };
+        float min_z = std::min(vertices[0].z(), std::min(vertices[1].z(), vertices[2].z()));
+        float max_z = std::max(vertices[0].z(), std::max(vertices[1].z(), vertices[2].z()));
+        
+        // intersect facet with cutting plane
+        IntersectionLine line;
+        int              idx_vertex_lowest = (vertices[1].z() == min_z) ? 1 : ((vertices[2].z() == min_z) ? 2 : 0);
+        FacetSliceType   slice_type        = slice_facet(scaled_z, vertices_scaled, m_its->indices[facet_idx], this->facets_edges.data() + facet_idx * 3, idx_vertex_lowest, min_z == max_z, &line);
+        if (slice_type != FacetSliceType::NoSlice) {
+            // Save intersection lines for generating correct triangulations.
+            if (line.edge_type == IntersectionLine::FacetEdgeType::Top) {
+                lower_lines.emplace_back(line);
+                lower_slice_vertices.emplace_back(line.a_id);
+                lower_slice_vertices.emplace_back(line.b_id);
+            } else if (line.edge_type == IntersectionLine::FacetEdgeType::Bottom) {
+                upper_lines.emplace_back(line);
+                upper_slice_vertices.emplace_back(line.a_id);
+                upper_slice_vertices.emplace_back(line.b_id);
+            } else if (line.edge_type == IntersectionLine::FacetEdgeType::General) {
+                lower_lines.emplace_back(line);
+                upper_lines.emplace_back(line);
+            }
+        }
+        
+        if (min_z > z || (min_z == z && max_z > z)) {
+            // facet is above the cut plane and does not belong to it
+            if (upper != nullptr)
+                upper->indices.emplace_back(facet);
+        } else if (max_z < z || (max_z == z && min_z < z)) {
+            // facet is below the cut plane and does not belong to it
+            if (lower != nullptr)
+                lower->indices.emplace_back(facet);
+        } else if (min_z < z && max_z > z) {
+            // Facet is cut by the slicing plane.
+            assert(slice_type == FacetSliceType::Slicing);
+            assert(line.edge_type == IntersectionLine::FacetEdgeType::General);
+
+            // look for the vertex on whose side of the slicing plane there are no other vertices
+            int isolated_vertex = 
+                (vertices[0].z() > z) == (vertices[1].z() > z) ? 2 :
+                (vertices[1].z() > z) == (vertices[2].z() > z) ? 0 : 1;
+            
+            // get vertices starting from the isolated one
+            int iv = isolated_vertex;
+            const stl_vertex &v0  = vertices[iv];
+            const int         iv0 = facet[iv];
+            if (++ iv == 3)
+                iv = 0;
+            const stl_vertex &v1  = vertices[iv];
+            const int         iv1 = facet[iv];
+            if (++ iv == 3)
+                iv = 0;
+            const stl_vertex &v2  = vertices[iv];
+            const int         iv2 = facet[iv];
+
+            // intersect v0-v1 and v2-v0 with cutting plane and make new vertices
+            auto new_vertex = [upper, lower, &upper_slice_vertices, &lower_slice_vertices](const Vec3f &a, const int ia, const Vec3f &b, const int ib, float z, float t) {
+                int iupper, ilower;
+                if (t <= 0.f)
+                    iupper = ilower = ia;
+                else if (t >= 1.f)
+                    iupper = ilower = ib;
+                else {
+                    const stl_vertex c = Vec3f(lerp(a.x(), b.x(), t), lerp(a.y(), b.y(), t), z);
+                    if (c == a)
+                        iupper = ilower = ia;
+                    else if (c == b)
+                        iupper = ilower = ib;
+                    else {
+                        // Insert a new vertex into upper / lower.
+                        if (upper) {
+                            iupper = int(upper->vertices.size());
+                            upper->vertices.emplace_back(c);
+                            upper_slice_vertices.emplace_back(iupper);
+                        }
+                        if (lower) {
+                            ilower = int(lower->vertices.size());
+                            lower->vertices.emplace_back(c);
+                            lower_slice_vertices.emplace_back(ilower);
+                        }
+                    }
+                }
+                return std::make_pair(iupper, ilower);
+            };
+            auto [iv0v1_upper, iv0v1_lower] = new_vertex(v1, iv1, v0, iv0, z, (z - v1.z()) / (v0.z() - v1.z()));
+            auto [iv2v0_upper, iv2v0_lower] = new_vertex(v2, iv2, v0, iv0, z, (z - v2.z()) / (v0.z() - v2.z()));
+            
+            if (v0(2) > z) {
+                if (upper != nullptr) 
+                    upper->indices.emplace_back(iv0, iv0v1_upper, iv2v0_upper);
+                if (lower != nullptr) {
+                    lower->indices.emplace_back(iv1, iv2, iv0v1_lower);
+                    lower->indices.emplace_back(iv2, iv2v0_lower, iv0v1_lower);
+                }
+            } else {
+                if (upper != nullptr) {
+                    upper->indices.emplace_back(iv1, iv2, iv0v1_upper);
+                    upper->indices.emplace_back(iv2, iv2v0_upper, iv0v1_upper);
+                }
+                if (lower != nullptr) 
+                    lower->indices.emplace_back(iv0, iv0v1_lower, iv2v0_lower);
+            }
+        }
+    }
+    
+    if (upper != nullptr)
+        triangulate_slice(*upper, upper_lines, upper_slice_vertices, z);
+
+    if (lower != nullptr)
+        triangulate_slice(*lower, lower_lines, lower_slice_vertices, z);
+}
+
+void TriangleMeshSlicer::cut(float z, TriangleMesh *upper_mesh, TriangleMesh *lower_mesh) const
+{
+    indexed_triangle_set upper, lower;
+    this->cut(z, upper_mesh ? &upper : nullptr, lower_mesh ? &lower : nullptr);
+    if (upper_mesh)
+        *upper_mesh = TriangleMesh(upper);
+    if (lower_mesh)
+        *lower_mesh = TriangleMesh(lower);
+}
+
+}
diff --git a/src/libslic3r/TriangleMeshSlicer.hpp b/src/libslic3r/TriangleMeshSlicer.hpp
new file mode 100644
index 000000000..258ffa04c
--- /dev/null
+++ b/src/libslic3r/TriangleMeshSlicer.hpp
@@ -0,0 +1,141 @@
+#ifndef slic3r_TriangleMeshSlicer_hpp_
+#define slic3r_TriangleMeshSlicer_hpp_
+
+#include "libslic3r.h"
+#include <admesh/stl.h>
+#include <functional>
+#include <vector>
+#include <boost/thread.hpp>
+#include "BoundingBox.hpp"
+#include "Line.hpp"
+#include "Point.hpp"
+#include "Polygon.hpp"
+#include "ExPolygon.hpp"
+
+namespace Slic3r {
+
+class TriangleMesh;
+
+enum class SlicingMode : uint32_t {
+    // Regular slicing, maintain all contours and their orientation.
+    Regular,
+    // Maintain all contours, orient all contours CCW, therefore all holes are being closed.
+    Positive,
+    // Orient all contours CCW and keep only the contour with the largest area.
+    // This mode is useful for slicing complex objects in vase mode.
+    PositiveLargestContour,
+};
+
+struct MeshSlicingParams
+{
+    SlicingMode   mode { SlicingMode::Regular };
+    // For vase mode: below this layer a different slicing mode will be used to produce a single contour.
+    // 0 = ignore.
+    size_t        slicing_mode_normal_below_layer { 0 };
+    // Mode to apply below slicing_mode_normal_below_layer. Ignored if slicing_mode_nromal_below_layer == 0.
+    SlicingMode   mode_below { SlicingMode::Regular };
+};
+
+struct MeshSlicingParamsExtended : public MeshSlicingParams
+{
+    // Morphological closing operation when creating output expolygons.
+    float         closing_radius { 0 };
+    // Positive offset applied when creating output expolygons.
+    float         extra_offset { 0 };
+    // Resolution for contour simplification.
+    // 0 = don't simplify.
+    double        resolution { 0 };
+    // Transformation of the object owning the ModelVolume.
+//    Transform3d         object_trafo;
+};
+
+class TriangleMeshSlicer
+{
+public:
+    using throw_on_cancel_callback_type = std::function<void()>;
+    TriangleMeshSlicer() = default;
+    TriangleMeshSlicer(const TriangleMesh *mesh) { this->init(mesh, []{}); }
+    TriangleMeshSlicer(const indexed_triangle_set *its) { this->init(its, []{}); }
+    void init(const TriangleMesh *mesh, throw_on_cancel_callback_type throw_on_cancel);
+    void init(const indexed_triangle_set *its, throw_on_cancel_callback_type);
+
+    void slice(
+        const std::vector<float>         &z,
+        const MeshSlicingParams          &params,
+        std::vector<Polygons>           *layers,
+        throw_on_cancel_callback_type    throw_on_cancel = []{}) const;
+
+    void slice(
+        // Where to slice.
+        const std::vector<float>         &z,
+        const MeshSlicingParamsExtended  &params,
+        std::vector<ExPolygons>         *layers,
+        throw_on_cancel_callback_type    throw_on_cancel = []{}) const;
+
+    void cut(float z, indexed_triangle_set *upper, indexed_triangle_set *lower) const;
+    void cut(float z, TriangleMesh* upper, TriangleMesh* lower) const;
+
+    void set_up_direction(const Vec3f& up);
+
+private:
+    const indexed_triangle_set *m_its { nullptr };
+//    const TriangleMesh      *mesh { nullptr };
+    // Map from a facet to an edge index.
+    std::vector<int>         facets_edges;
+    // Scaled copy of this->mesh->stl.v_shared
+    std::vector<stl_vertex>  v_scaled_shared;
+    // Quaternion that will be used to rotate every facet before the slicing
+    Eigen::Quaternion<float, Eigen::DontAlign> m_quaternion;
+    // Whether or not the above quaterion should be used
+    bool                     m_use_quaternion = false;
+};
+
+inline void slice_mesh(
+    const TriangleMesh                               &mesh,
+    const std::vector<float>                         &z,
+    std::vector<Polygons>                            &layers,
+    TriangleMeshSlicer::throw_on_cancel_callback_type thr = []{})
+{
+    if (! mesh.empty()) {
+        TriangleMeshSlicer slicer(&mesh);
+        slicer.slice(z, MeshSlicingParams{}, &layers, thr);
+    }
+}
+
+inline void slice_mesh(
+    const TriangleMesh                               &mesh,
+    const std::vector<float>                         &z,
+    const MeshSlicingParamsExtended                  &params,
+    std::vector<ExPolygons>                          &layers,
+    TriangleMeshSlicer::throw_on_cancel_callback_type thr = []{})
+{
+    if (! mesh.empty()) {
+        TriangleMeshSlicer slicer(&mesh);
+        slicer.slice(z, params, &layers, thr);
+    }
+}
+
+inline void slice_mesh(
+    const TriangleMesh                               &mesh,
+    const std::vector<float>                         &z,
+    float                                             closing_radius,
+    std::vector<ExPolygons>                          &layers,
+    TriangleMeshSlicer::throw_on_cancel_callback_type thr = []{})
+{
+    MeshSlicingParamsExtended params;
+    params.closing_radius = closing_radius;
+    slice_mesh(mesh, z, params, layers);
+}
+
+inline void slice_mesh(
+    const TriangleMesh                               &mesh,
+    const std::vector<float>                         &z,
+    std::vector<ExPolygons>                          &layers,
+    TriangleMeshSlicer::throw_on_cancel_callback_type thr = []{})
+{
+    slice_mesh(mesh, z, MeshSlicingParamsExtended{}, layers);
+}
+
+}
+
+#endif // slic3r_TriangleMeshSlicer_hpp_
diff --git a/src/slic3r/GUI/BackgroundSlicingProcess.hpp b/src/slic3r/GUI/BackgroundSlicingProcess.hpp
index dc24ec0ad..12bf6fe02 100644
--- a/src/slic3r/GUI/BackgroundSlicingProcess.hpp
+++ b/src/slic3r/GUI/BackgroundSlicingProcess.hpp
@@ -5,6 +5,8 @@
 #include <condition_variable>
 #include <mutex>
 
+#include <boost/thread.hpp>
+
 #include <wx/event.h>
 
 #include "libslic3r/PrintBase.hpp"
diff --git a/src/slic3r/GUI/DoubleSlider.cpp b/src/slic3r/GUI/DoubleSlider.cpp
index a10acb2e3..0c30b0cb4 100644
--- a/src/slic3r/GUI/DoubleSlider.cpp
+++ b/src/slic3r/GUI/DoubleSlider.cpp
@@ -254,7 +254,7 @@ void Control::SetMaxValue(const int max_value)
 void Control::SetSliderValues(const std::vector<double>& values)
 {
     m_values = values;
-    m_ruler.count = std::count(m_values.begin(), m_values.end(), m_values.front());
+    m_ruler.init(m_values);
 }
 
 void Control::draw_scroll_line(wxDC& dc, const int lower_pos, const int higher_pos)
@@ -1023,8 +1023,23 @@ void Control::draw_colored_band(wxDC& dc)
     }
 }
 
+void Control::Ruler::init(const std::vector<double>& values)
+{
+    max_values.clear();
+    max_values.reserve(std::count(values.begin(), values.end(), values.front()));
+
+    auto it = std::find(values.begin() + 1, values.end(), values.front());
+    while (it != values.end()) {
+        max_values.push_back(*(it - 1));
+        it = std::find(it + 1, values.end(), values.front());
+    }
+    max_values.push_back(*(it - 1));
+}
+
 void Control::Ruler::update(wxWindow* win, const std::vector<double>& values, double scroll_step)
 {
+    if (values.empty())
+        return;
     int DPI = GUI::get_dpi_for_window(win);
     int pixels_per_sm = lround((double)(DPI) * 5.0/25.4);
 
@@ -1035,7 +1050,7 @@ void Control::Ruler::update(wxWindow* win, const std::vector<double>& values, do
 
     int pow = -2;
     int step = 0;
-    auto end_it = count == 1 ? values.end() : values.begin() + lround(values.size() / count);
+    auto end_it = std::find(values.begin() + 1, values.end(), values.front());
 
     while (pow < 3) {
         for (int istep : {1, 2, 5}) {
@@ -1099,7 +1114,7 @@ void Control::draw_ruler(wxDC& dc)
         double short_tick = std::nan("");
         int tick = 0;
         double value = 0.0;
-        int sequence = 0;
+        size_t sequence = 0;
 
         int prev_y_pos = -1;
         wxCoord label_height = dc.GetMultiLineTextExtent("0").y - 2;
@@ -1107,7 +1122,7 @@ void Control::draw_ruler(wxDC& dc)
 
         while (tick <= m_max_value) {
             value += m_ruler.long_step;
-            if (value > m_values.back() && sequence < m_ruler.count) {
+            if (value > m_ruler.max_values[sequence] && sequence < m_ruler.count()) {
                 value = m_ruler.long_step;
                 for (; tick < values_size; tick++)
                     if (m_values[tick] < value)
@@ -1140,7 +1155,7 @@ void Control::draw_ruler(wxDC& dc)
 
             draw_short_ticks(dc, short_tick, tick);
 
-            if (value == m_values.back() && sequence < m_ruler.count) {
+            if (value == m_ruler.max_values[sequence] && sequence < m_ruler.count()) {
                 value = 0.0;
                 sequence++;
                 tick++;
diff --git a/src/slic3r/GUI/DoubleSlider.hpp b/src/slic3r/GUI/DoubleSlider.hpp
index 400265885..c1766f83d 100644
--- a/src/slic3r/GUI/DoubleSlider.hpp
+++ b/src/slic3r/GUI/DoubleSlider.hpp
@@ -424,10 +424,13 @@ private:
     struct Ruler {
         double long_step;
         double short_step;
-        int count { 1 }; // > 1 for sequential print
+        std::vector<double> max_values;// max value for each object/instance in sequence print
+                                       // > 1 for sequential print
 
+        void init(const std::vector<double>& values);
         void update(wxWindow* win, const std::vector<double>& values, double scroll_step);
         bool is_ok() { return long_step > 0 && short_step > 0; }
+        size_t count() { return max_values.size(); }
     } m_ruler;
 };
 
diff --git a/src/slic3r/GUI/MeshUtils.cpp b/src/slic3r/GUI/MeshUtils.cpp
index f9ccfd0d6..2f575b504 100644
--- a/src/slic3r/GUI/MeshUtils.cpp
+++ b/src/slic3r/GUI/MeshUtils.cpp
@@ -2,6 +2,7 @@
 
 #include "libslic3r/Tesselate.hpp"
 #include "libslic3r/TriangleMesh.hpp"
+#include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/ClipperUtils.hpp"
 
 #include "slic3r/GUI/Camera.hpp"
@@ -83,16 +84,17 @@ void MeshClipper::recalculate_triangles()
     // Now do the cutting
     std::vector<ExPolygons> list_of_expolys;
     m_tms->set_up_direction(up.cast<float>());
-    m_tms->slice(std::vector<float>{height_mesh}, SlicingMode::Regular, 0.f, &list_of_expolys, [](){});
+    m_tms->slice(std::vector<float>{height_mesh}, MeshSlicingParamsExtended{}, &list_of_expolys);
 
     if (m_negative_mesh && !m_negative_mesh->empty()) {
         TriangleMeshSlicer negative_tms{m_negative_mesh};
         negative_tms.set_up_direction(up.cast<float>());
 
         std::vector<ExPolygons> neg_polys;
-        negative_tms.slice(std::vector<float>{height_mesh}, SlicingMode::Regular, 0.f, &neg_polys, [](){});
+        negative_tms.slice(std::vector<float>{height_mesh}, MeshSlicingParamsExtended{}, &neg_polys);
         list_of_expolys.front() = diff_ex(list_of_expolys.front(), neg_polys.front());
     }
+   
     m_triangles2d = triangulate_expolygons_2f(list_of_expolys[0], m_trafo.get_matrix().matrix().determinant() < 0.);
 
     // Rotate the cut into world coords:
diff --git a/src/slic3r/GUI/MeshUtils.hpp b/src/slic3r/GUI/MeshUtils.hpp
index 09caf199b..8a430e322 100644
--- a/src/slic3r/GUI/MeshUtils.hpp
+++ b/src/slic3r/GUI/MeshUtils.hpp
@@ -3,6 +3,7 @@
 
 #include "libslic3r/Point.hpp"
 #include "libslic3r/Geometry.hpp"
+#include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/SLA/IndexedMesh.hpp"
 #include "admesh/stl.h"
 
@@ -12,9 +13,6 @@
 
 namespace Slic3r {
 
-class TriangleMesh;
-class TriangleMeshSlicer;
-
 namespace GUI {
 
 struct Camera;
diff --git a/src/slic3r/Utils/FixModelByWin10.cpp b/src/slic3r/Utils/FixModelByWin10.cpp
index 7933a1d69..dead35267 100644
--- a/src/slic3r/Utils/FixModelByWin10.cpp
+++ b/src/slic3r/Utils/FixModelByWin10.cpp
@@ -27,6 +27,7 @@
 #include <boost/filesystem.hpp>
 #include <boost/nowide/convert.hpp>
 #include <boost/nowide/cstdio.hpp>
+#include <boost/thread.hpp>
 
 #include "libslic3r/Model.hpp"
 #include "libslic3r/Print.hpp"
diff --git a/tests/fff_print/test_trianglemesh.cpp b/tests/fff_print/test_trianglemesh.cpp
index 233b0e515..a066215c2 100644
--- a/tests/fff_print/test_trianglemesh.cpp
+++ b/tests/fff_print/test_trianglemesh.cpp
@@ -1,6 +1,7 @@
 #include <catch2/catch.hpp>
 
 #include "libslic3r/TriangleMesh.hpp"
+#include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/Point.hpp"
 #include "libslic3r/Config.hpp"
 #include "libslic3r/Model.hpp"
diff --git a/tests/libslic3r/test_marchingsquares.cpp b/tests/libslic3r/test_marchingsquares.cpp
index 1a4b1fb72..66a0060fe 100644
--- a/tests/libslic3r/test_marchingsquares.cpp
+++ b/tests/libslic3r/test_marchingsquares.cpp
@@ -13,6 +13,7 @@
 #include <libslic3r/SVG.hpp>
 #include <libslic3r/ClipperUtils.hpp>
 
+#include <libslic3r/TriangleMeshSlicer.hpp>
 #include <libslic3r/TriangulateWall.hpp>
 #include <libslic3r/Tesselate.hpp>
 #include <libslic3r/SlicesToTriangleMesh.hpp>
@@ -320,7 +321,7 @@ static void recreate_object_from_rasters(const std::string &objname, float lh) {
     bb = mesh.bounding_box();
     
     std::vector<ExPolygons> layers;
-    slice_mesh(mesh, grid(float(bb.min.z()) + lh, float(bb.max.z()), lh), layers, 0.f, []{});
+    slice_mesh(mesh, grid(float(bb.min.z()) + lh, float(bb.max.z()), lh), layers);
     
     sla::RasterBase::Resolution res{2560, 1440};
     double                      disp_w = 120.96;
diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp
index 1f98463cc..75748dd34 100644
--- a/tests/sla_print/sla_print_tests.cpp
+++ b/tests/sla_print/sla_print_tests.cpp
@@ -6,6 +6,7 @@
 
 #include "sla_test_utils.hpp"
 
+#include <libslic3r/TriangleMeshSlicer.hpp>
 #include <libslic3r/SLA/SupportTreeMesher.hpp>
 #include <libslic3r/SLA/Concurrency.hpp>
 
@@ -48,9 +49,7 @@ TEST_CASE("Support point generator should be deterministic if seeded",
     sla::SupportPointGenerator::Config autogencfg;
     autogencfg.head_diameter = float(2 * supportcfg.head_front_radius_mm);
     sla::SupportPointGenerator point_gen{emesh, autogencfg, [] {}, [](int) {}};
-    
-    TriangleMeshSlicer slicer{&mesh};
-    
+        
     auto   bb      = mesh.bounding_box();
     double zmin    = bb.min.z();
     double zmax    = bb.max.z();
@@ -59,7 +58,7 @@ TEST_CASE("Support point generator should be deterministic if seeded",
     
     auto slicegrid = grid(float(gnd), float(zmax), layer_h);
     std::vector<ExPolygons> slices;
-    slicer.slice(slicegrid, SlicingMode::Regular, CLOSING_RADIUS, &slices, []{});
+    slice_mesh(mesh, slicegrid, CLOSING_RADIUS, slices);
     
     point_gen.seed(0);
     point_gen.execute(slices, slicegrid);
diff --git a/tests/sla_print/sla_test_utils.cpp b/tests/sla_print/sla_test_utils.cpp
index 1ec890beb..fdf77466b 100644
--- a/tests/sla_print/sla_test_utils.cpp
+++ b/tests/sla_print/sla_test_utils.cpp
@@ -1,4 +1,5 @@
 #include "sla_test_utils.hpp"
+#include "libslic3r/TriangleMeshSlicer.hpp"
 #include "libslic3r/SLA/AGGRaster.hpp"
 
 void test_support_model_collision(const std::string          &obj_filename,
@@ -94,8 +95,6 @@ void test_supports(const std::string          &obj_filename,
         mesh.require_shared_vertices();
     }
     
-    TriangleMeshSlicer slicer{&mesh};
-    
     auto   bb      = mesh.bounding_box();
     double zmin    = bb.min.z();
     double zmax    = bb.max.z();
@@ -103,7 +102,7 @@ void test_supports(const std::string          &obj_filename,
     auto   layer_h = 0.05f;
     
     out.slicegrid = grid(float(gnd), float(zmax), layer_h);
-    slicer.slice(out.slicegrid, SlicingMode::Regular, CLOSING_RADIUS, &out.model_slices, []{});
+    slice_mesh(mesh, out.slicegrid, CLOSING_RADIUS, out.model_slices);
     sla::cut_drainholes(out.model_slices, out.slicegrid, CLOSING_RADIUS, drainholes, []{});
     
     // Create the special index-triangle mesh with spatial indexing which
@@ -470,7 +469,7 @@ sla::SupportPoints calc_support_pts(
     std::vector<ExPolygons> slices;
     auto                    bb      = cast<float>(mesh.bounding_box());
     std::vector<float>      heights = grid(bb.min.z(), bb.max.z(), 0.1f);
-    slice_mesh(mesh, heights, slices, CLOSING_RADIUS, [] {});
+    slice_mesh(mesh, heights, CLOSING_RADIUS, slices);
 
     // Prepare the support point calculator
     sla::IndexedMesh emesh{mesh};
diff --git a/xs/xsp/TriangleMesh.xsp b/xs/xsp/TriangleMesh.xsp
index 230f8b2a5..377bf7b5e 100644
--- a/xs/xsp/TriangleMesh.xsp
+++ b/xs/xsp/TriangleMesh.xsp
@@ -3,6 +3,7 @@
 %{
 #include <xsinit.h>
 #include "libslic3r/TriangleMesh.hpp"
+#include "libslic3r/TriangleMeshSlicer.hpp"
 %}
 
 %name{Slic3r::TriangleMesh} class TriangleMesh {
@@ -181,8 +182,7 @@ TriangleMesh::slice(z)
         std::vector<float> z_f = cast<float>(z);
         
         std::vector<ExPolygons> layers;
-        TriangleMeshSlicer mslicer(THIS);
-        mslicer.slice(z_f, SlicingMode::Regular, 0.049f, &layers, [](){});
+        slice_mesh(*THIS, z_f, 0.049f, layers);
         
         AV* layers_av = newAV();
         size_t len = layers.size();