TriangleMeshSlicer: Got rid of admesh!

11 files changed, 408 insertions, 410 deletions

diff --git a/src/libslic3r/Model.cpp b/src/libslic3r/Model.cpp
index 4ddcbac39..55eb6b995 100644
--- a/src/libslic3r/Model.cpp
+++ b/src/libslic3r/Model.cpp
@@ -1213,32 +1213,32 @@ ModelObjectPtrs ModelObject::cut(size_t instance, coordf_t z, bool keep_upper, b
         }
         else if (! volume->mesh().empty()) {
             
-            TriangleMesh upper_mesh, lower_mesh;
-
             // Transform the mesh by the combined transformation matrix.
             // Flip the triangles in case the composite transformation is left handed.
 			TriangleMesh mesh(volume->mesh());
 			mesh.transform(instance_matrix * volume_matrix, true);
 			volume->reset_mesh();
-            
-            mesh.require_shared_vertices();
-            
-            // Perform cut
-            TriangleMeshSlicer tms(&mesh);
-            tms.cut(float(z), &upper_mesh, &lower_mesh);
-
             // Reset volume transformation except for offset
             const Vec3d offset = volume->get_offset();
             volume->set_transformation(Geometry::Transformation());
             volume->set_offset(offset);
 
-            if (keep_upper) {
-                upper_mesh.repair();
-                upper_mesh.reset_repair_stats();
-            }
-            if (keep_lower) {
-                lower_mesh.repair();
-                lower_mesh.reset_repair_stats();
+            // Perform cut
+            TriangleMesh upper_mesh, lower_mesh;
+            {
+                indexed_triangle_set upper_its, lower_its;
+                mesh.require_shared_vertices();
+                cut_mesh(mesh.its, float(z), &upper_its, &lower_its);
+                if (keep_upper) {
+                    upper_mesh = TriangleMesh(upper_its);
+                    upper_mesh.repair();
+                    upper_mesh.reset_repair_stats();
+                }
+                if (keep_lower) {
+                    lower_mesh = TriangleMesh(lower_its);
+                    lower_mesh.repair();
+                    lower_mesh.reset_repair_stats();
+                }
             }
 
             if (keep_upper && upper_mesh.facets_count() > 0) {
diff --git a/src/libslic3r/Print.hpp b/src/libslic3r/Print.hpp
index 3fdc49db8..edd8e19f1 100644
--- a/src/libslic3r/Print.hpp
+++ b/src/libslic3r/Print.hpp
@@ -8,6 +8,7 @@
 #include "Flow.hpp"
 #include "Point.hpp"
 #include "Slicing.hpp"
+#include "TriangleMeshSlicer.hpp"
 #include "GCode/ToolOrdering.hpp"
 #include "GCode/WipeTower.hpp"
 #include "GCode/ThumbnailData.hpp"
@@ -24,7 +25,6 @@ class Print;
 class PrintObject;
 class ModelObject;
 class GCode;
-enum class SlicingMode : uint32_t;
 class Layer;
 class SupportLayer;
 
@@ -345,18 +345,18 @@ private:
     // so that next call to make_perimeters() performs a union() before computing loops
     bool                    				m_typed_slices = false;
 
-    std::vector<ExPolygons> slice_region(size_t region_id, const std::vector<float> &z, SlicingMode mode, size_t slicing_mode_normal_below_layer, SlicingMode mode_below) const;
-    std::vector<ExPolygons> slice_region(size_t region_id, const std::vector<float> &z, SlicingMode mode) const
+    std::vector<ExPolygons> slice_region(size_t region_id, const std::vector<float> &z, MeshSlicingParams::SlicingMode mode, size_t slicing_mode_normal_below_layer, MeshSlicingParams::SlicingMode mode_below) const;
+    std::vector<ExPolygons> slice_region(size_t region_id, const std::vector<float> &z, MeshSlicingParams::SlicingMode mode) const
         { return this->slice_region(region_id, z, mode, 0, mode); }
     std::vector<ExPolygons> slice_modifiers(size_t region_id, const std::vector<float> &z) const;
     std::vector<ExPolygons> slice_volumes(
         const std::vector<float> &z, 
-        SlicingMode mode, size_t slicing_mode_normal_below_layer, SlicingMode mode_below, 
+        MeshSlicingParams::SlicingMode mode, size_t slicing_mode_normal_below_layer, MeshSlicingParams::SlicingMode mode_below,
         const std::vector<const ModelVolume*> &volumes) const;
-    std::vector<ExPolygons> slice_volumes(const std::vector<float> &z, SlicingMode mode, const std::vector<const ModelVolume*> &volumes) const
+    std::vector<ExPolygons> slice_volumes(const std::vector<float> &z, MeshSlicingParams::SlicingMode mode, const std::vector<const ModelVolume*> &volumes) const
         { return this->slice_volumes(z, mode, 0, mode, volumes); }
-    std::vector<ExPolygons> slice_volume(const std::vector<float> &z, SlicingMode mode, const ModelVolume &volume) const;
-    std::vector<ExPolygons> slice_volume(const std::vector<float> &z, const std::vector<t_layer_height_range> &ranges, SlicingMode mode, const ModelVolume &volume) const;
+    std::vector<ExPolygons> slice_volume(const std::vector<float> &z, MeshSlicingParams::SlicingMode mode, const ModelVolume &volume) const;
+    std::vector<ExPolygons> slice_volume(const std::vector<float> &z, const std::vector<t_layer_height_range> &ranges, MeshSlicingParams::SlicingMode mode, const ModelVolume &volume) const;
 };
 
 struct WipeTowerData
diff --git a/src/libslic3r/PrintObject.cpp b/src/libslic3r/PrintObject.cpp
index b77b13345..20c45f3cb 100644
--- a/src/libslic3r/PrintObject.cpp
+++ b/src/libslic3r/PrintObject.cpp
@@ -1799,7 +1799,7 @@ void PrintObject::_slice(const std::vector<coordf_t> &layer_height_profile)
     bool clipped  = false;
     bool upscaled = false;
     bool spiral_vase  = this->print()->config().spiral_vase;
-    auto slicing_mode = spiral_vase ? SlicingMode::PositiveLargestContour : SlicingMode::Regular;
+    auto slicing_mode = spiral_vase ? MeshSlicingParams::SlicingMode::PositiveLargestContour : MeshSlicingParams::SlicingMode::Regular;
     if (! has_z_ranges && (! m_config.clip_multipart_objects.value || all_volumes_single_region >= 0)) {
         // Cheap path: Slice regions without mutual clipping.
         // The cheap path is possible if no clipping is allowed or if slicing volumes of just a single region.
@@ -1815,7 +1815,7 @@ void PrintObject::_slice(const std::vector<coordf_t> &layer_height_profile)
                 for (; slicing_mode_normal_below_layer < slice_zs.size() && slice_zs[slicing_mode_normal_below_layer] < config.bottom_solid_min_thickness - EPSILON;
                     ++ slicing_mode_normal_below_layer);
             }
-            std::vector<ExPolygons> expolygons_by_layer = this->slice_region(region_id, slice_zs, slicing_mode, slicing_mode_normal_below_layer, SlicingMode::Regular);
+            std::vector<ExPolygons> expolygons_by_layer = this->slice_region(region_id, slice_zs, slicing_mode, slicing_mode_normal_below_layer, MeshSlicingParams::SlicingMode::Regular);
             m_print->throw_if_canceled();
             BOOST_LOG_TRIVIAL(debug) << "Slicing objects - append slices " << region_id << " start";
             for (size_t layer_id = 0; layer_id < expolygons_by_layer.size(); ++ layer_id)
@@ -2060,7 +2060,7 @@ end:
 }
 
 // To be used only if there are no layer span specific configurations applied, which would lead to z ranges being generated for this region.
-std::vector<ExPolygons> PrintObject::slice_region(size_t region_id, const std::vector<float> &z, SlicingMode mode, size_t slicing_mode_normal_below_layer, SlicingMode mode_below) const
+std::vector<ExPolygons> PrintObject::slice_region(size_t region_id, const std::vector<float> &z, MeshSlicingParams::SlicingMode mode, size_t slicing_mode_normal_below_layer, MeshSlicingParams::SlicingMode mode_below) const
 {
 	std::vector<const ModelVolume*> volumes;
     if (region_id < m_region_volumes.size()) {
@@ -2120,7 +2120,7 @@ std::vector<ExPolygons> PrintObject::slice_modifiers(size_t region_id, const std
 					if (volume->is_modifier())
 						volumes.emplace_back(volume);
 				}
-				out = this->slice_volumes(slice_zs, SlicingMode::Regular, volumes);
+				out = this->slice_volumes(slice_zs, MeshSlicingParams::SlicingMode::Regular, volumes);
 			} else {
 				// Some modifier in this region was split to layer spans.
 				std::vector<char> merge;
@@ -2138,7 +2138,7 @@ std::vector<ExPolygons> PrintObject::slice_modifiers(size_t region_id, const std
 							for (; j < volumes_and_ranges.volumes.size() && volume_id == volumes_and_ranges.volumes[j].volume_idx; ++ j)
 								ranges.emplace_back(volumes_and_ranges.volumes[j].layer_height_range);
 			                // slicing in parallel
-			                std::vector<ExPolygons> this_slices = this->slice_volume(slice_zs, ranges, SlicingMode::Regular, *model_volume);
+			                std::vector<ExPolygons> this_slices = this->slice_volume(slice_zs, ranges, MeshSlicingParams::SlicingMode::Regular, *model_volume);
                             // Variable this_slices could be empty if no value of slice_zs is within any of the ranges of this volume.
 			                if (out.empty()) {
 			                	out = std::move(this_slices);
@@ -2179,7 +2179,7 @@ std::vector<ExPolygons> PrintObject::slice_support_volumes(const ModelVolumeType
     zs.reserve(this->layers().size());
     for (const Layer *l : this->layers())
         zs.emplace_back((float)l->slice_z);
-    return this->slice_volumes(zs, SlicingMode::Regular, volumes);
+    return this->slice_volumes(zs, MeshSlicingParams::SlicingMode::Regular, volumes);
 }
 
 //FIXME The admesh repair function may break the face connectivity, rather refresh it here as the slicing code relies on it.
@@ -2194,7 +2194,7 @@ static void fix_mesh_connectivity(TriangleMesh &mesh)
 
 std::vector<ExPolygons> PrintObject::slice_volumes(
     const std::vector<float> &z, 
-    SlicingMode mode, size_t slicing_mode_normal_below_layer, SlicingMode mode_below, 
+    MeshSlicingParams::SlicingMode mode, size_t slicing_mode_normal_below_layer, MeshSlicingParams::SlicingMode mode_below,
     const std::vector<const ModelVolume*> &volumes) const
 {
     std::vector<ExPolygons> layers;
@@ -2219,18 +2219,17 @@ std::vector<ExPolygons> PrintObject::slice_volumes(
             mesh.translate(- unscale<float>(m_center_offset.x()), - unscale<float>(m_center_offset.y()), 0);
             // perform actual slicing
             const Print *print = this->print();
-            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();
-            MeshSlicingParamsExtended params { { mode, slicing_mode_normal_below_layer, mode_below }, float(m_config.slice_closing_radius.value) };
-			slice_mesh(mesh, z, params, layers, callback);
+            MeshSlicingParamsEx params { { mode, slicing_mode_normal_below_layer, mode_below }, float(m_config.slice_closing_radius.value) };
+            layers = slice_mesh_ex(mesh.its, z, params, [print]() { print->throw_if_canceled(); });
             m_print->throw_if_canceled();
         }
     }
     return layers;
 }
 
-std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, SlicingMode mode, const ModelVolume &volume) const
+std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, MeshSlicingParams::SlicingMode mode, const ModelVolume &volume) const
 {
     std::vector<ExPolygons> layers;
     if (! z.empty()) {
@@ -2246,13 +2245,12 @@ std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, S
 	        mesh.translate(- unscale<float>(m_center_offset.x()), - unscale<float>(m_center_offset.y()), 0);
 	        // perform actual slicing
 	        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();
-            MeshSlicingParamsExtended params;
+            MeshSlicingParamsEx params;
             params.mode = mode;
             params.closing_radius = float(m_config.slice_closing_radius.value);
-	        slice_mesh(mesh, z, params, layers, callback);
+            layers = slice_mesh_ex(mesh.its, z, params, [print](){ print->throw_if_canceled(); });
 	        m_print->throw_if_canceled();
 	    }
 	}
@@ -2260,7 +2258,7 @@ std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, S
 }
 
 // Filter the zs not inside the ranges. The ranges are closed at the bottom and open at the top, they are sorted lexicographically and non overlapping.
-std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, const std::vector<t_layer_height_range> &ranges, SlicingMode mode, const ModelVolume &volume) const
+std::vector<ExPolygons> PrintObject::slice_volume(const std::vector<float> &z, const std::vector<t_layer_height_range> &ranges, MeshSlicingParams::SlicingMode mode, const ModelVolume &volume) const
 {
 	std::vector<ExPolygons> out;
 	if (! z.empty() && ! ranges.empty()) {
diff --git a/src/libslic3r/SLA/Hollowing.cpp b/src/libslic3r/SLA/Hollowing.cpp
index f55a10178..1d3016bde 100644
--- a/src/libslic3r/SLA/Hollowing.cpp
+++ b/src/libslic3r/SLA/Hollowing.cpp
@@ -296,9 +296,8 @@ void cut_drainholes(std::vector<ExPolygons> & obj_slices,
     if (mesh.empty()) return;
     
     mesh.require_shared_vertices();
-    
-    std::vector<ExPolygons> hole_slices;
-    slice_mesh(mesh, slicegrid, closing_radius, hole_slices, thr);
+ 
+    std::vector<ExPolygons> hole_slices = slice_mesh_ex(mesh.its, slicegrid, closing_radius, 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 658ecf6d8..bf2b4cf54 100644
--- a/src/libslic3r/SLA/Pad.cpp
+++ b/src/libslic3r/SLA/Pad.cpp
@@ -478,8 +478,8 @@ void pad_blueprint(const TriangleMesh &      mesh,
 {
     if (mesh.empty()) return;
 
-    auto out = reserve_vector<ExPolygons>(heights.size());
-    slice_mesh(mesh, heights, out, thrfn);
+    assert(mesh.has_shared_vertices());
+    std::vector<ExPolygons> out = slice_mesh_ex(mesh.its, heights, thrfn);
 
     size_t count = 0;
     for(auto& o : out) count += o.size();
diff --git a/src/libslic3r/SLA/SupportTree.cpp b/src/libslic3r/SLA/SupportTree.cpp
index 45e90a704..14a4dc360 100644
--- a/src/libslic3r/SLA/SupportTree.cpp
+++ b/src/libslic3r/SLA/SupportTree.cpp
@@ -45,7 +45,8 @@ std::vector<ExPolygons> SupportTree::slice(
 
     if (!sup_mesh.empty()) {
         slices.emplace_back();
-        slice_mesh(sup_mesh, grid, cr, slices.back(), ctl().cancelfn);
+        assert(sup_mesh.has_shared_vertices());
+        slices.back() = slice_mesh_ex(sup_mesh.its, grid, cr, ctl().cancelfn);
     }
 
     if (!pad_mesh.empty()) {
@@ -58,7 +59,8 @@ 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));
 
-        slice_mesh(pad_mesh, padgrid, cr, slices.back(), ctl().cancelfn);
+        assert(pad_mesh.has_shared_vertices());
+        slices.back() = slice_mesh_ex(pad_mesh.its, padgrid, cr, ctl().cancelfn);
     }
 
     size_t len = grid.size();
diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp
index d2c3c06d2..46064c55b 100644
--- a/src/libslic3r/SLAPrintSteps.cpp
+++ b/src/libslic3r/SLAPrintSteps.cpp
@@ -475,7 +475,8 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po)
     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;
-    slice_mesh(mesh, slice_grid, closing_r, po.m_model_slices, thr);
+    assert(mesh.has_shared_vertices());
+    po.m_model_slices = slice_mesh_ex(mesh.its, slice_grid, closing_r, thr);
 
     sla::Interior *interior = po.m_hollowing_data ?
                                   po.m_hollowing_data->interior.get() :
@@ -485,8 +486,7 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po)
         TriangleMesh interiormesh = sla::get_mesh(*interior);
         interiormesh.repaired = false;
         interiormesh.repair(true);
-        std::vector<ExPolygons> interior_slices;
-        slice_mesh(interiormesh, slice_grid, closing_r, interior_slices, thr);
+        std::vector<ExPolygons> interior_slices = slice_mesh_ex(interiormesh.its, slice_grid, closing_r, thr);
 
         sla::ccr::for_each(size_t(0), interior_slices.size(),
                            [&po, &interior_slices] (size_t i) {
diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp
index 9ce2d3a7f..2c35c76a2 100644
--- a/src/libslic3r/TriangleMesh.cpp
+++ b/src/libslic3r/TriangleMesh.cpp
@@ -615,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());
-    std::vector<ExPolygons> layers;
-    slice_mesh(*this, z_f, 0.0004f, layers);
-    return layers;
+    assert(this->has_shared_vertices());
+    return slice_mesh_ex(this->its, z_f, 0.0004f);
 }
 
 void TriangleMesh::require_shared_vertices()
@@ -686,11 +685,12 @@ std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangl
     return index;
 }
 
-// Map from a facet edge to a neighbor face index or -1 if no neighbor exists.
+// Map from a face edge to a unique edge identifier or -1 if no neighbor exists.
+// Two neighbor faces share a unique edge identifier even if they are flipped.
 template<typename ThrowOnCancelCallback>
-static inline std::vector<int> create_face_neighbors_index_impl(const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
+static inline std::vector<Vec3i> create_face_neighbors_index_impl(const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
 {
-    std::vector<int> out(its.indices.size() * 3, -1);
+    std::vector<Vec3i> out(its.indices.size(), Vec3i(-1, -1, -1));
 
     // Create a mapping from triangle edge into face.
     struct EdgeToFace {
@@ -754,10 +754,10 @@ static inline std::vector<int> create_face_neighbors_index_impl(const indexed_tr
                 }
         }
         // Assign an edge index to the 1st face.
-        out[edge_i.face * 3 + std::abs(edge_i.face_edge) - 1] = num_edges;
+        out[edge_i.face](std::abs(edge_i.face_edge) - 1) = num_edges;
         if (found) {
             EdgeToFace &edge_j = edges_map[j];
-            out[edge_j.face * 3 + std::abs(edge_j.face_edge) - 1] = num_edges;
+            out[edge_j.face](std::abs(edge_j.face_edge) - 1) = num_edges;
             // Mark the edge as connected.
             edge_j.face = -1;
         }
@@ -769,12 +769,12 @@ static inline std::vector<int> create_face_neighbors_index_impl(const indexed_tr
     return out;
 }
 
-std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its)
+std::vector<Vec3i> create_face_neighbors_index(const indexed_triangle_set &its)
 {
     return create_face_neighbors_index_impl(its, [](){});
 }
 
-std::vector<int> create_face_neighbors_index(const indexed_triangle_set &its, std::function<void()> throw_on_cancel_callback)
+std::vector<Vec3i> create_face_neighbors_index(const indexed_triangle_set &its, std::function<void()> throw_on_cancel_callback)
 {
     return create_face_neighbors_index_impl(its, throw_on_cancel_callback);
 }
diff --git a/src/libslic3r/TriangleMesh.hpp b/src/libslic3r/TriangleMesh.hpp
index dd1a76156..2be822350 100644
--- a/src/libslic3r/TriangleMesh.hpp
+++ b/src/libslic3r/TriangleMesh.hpp
@@ -93,9 +93,11 @@ private:
 // vertex.
 std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangle_set &its);
 
-// 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);
+// Map from a face edge to a unique edge identifier or -1 if no neighbor exists.
+// Two neighbor faces share a unique edge identifier even if they are flipped.
+// Used for chaining slice lines into polygons.
+std::vector<Vec3i> create_face_neighbors_index(const indexed_triangle_set &its);
+std::vector<Vec3i> 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.
diff --git a/src/libslic3r/TriangleMeshSlicer.cpp b/src/libslic3r/TriangleMeshSlicer.cpp
index 657a59b62..06061b034 100644
--- a/src/libslic3r/TriangleMeshSlicer.cpp
+++ b/src/libslic3r/TriangleMeshSlicer.cpp
@@ -85,7 +85,7 @@ public:
 
     // feGeneral, feTop, feBottom, feHorizontal
     FacetEdgeType   edge_type { FacetEdgeType::General };
-    // Used by TriangleMeshSlicer::slice() to skip duplicate edges.
+    // Used to skip duplicate edges.
     enum {
         // Triangle edge added, because it has no neighbor.
         EDGE0_NO_NEIGHBOR   = 0x001,
@@ -112,8 +112,15 @@ enum class FacetSliceType {
 
 // 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)
+    // Z height of the slice in XY plane. Scaled or unscaled (same as vertices[].z()).
+    float                                slice_z,
+    // 3 vertices of the triangle, XY scaled. Z scaled or unscaled (same as slice_z).
+    const stl_vertex                    *vertices,
+    const stl_triangle_vertex_indices   &indices,
+    const Vec3i                         &edge_neighbor,
+    const int                            idx_vertex_lowest,
+    const bool                           horizontal,
+    IntersectionLine                    &line_out)
 {
     IntersectionPoint points[3];
     size_t            num_points = 0;
@@ -129,7 +136,7 @@ static FacetSliceType slice_facet(
         {
             int   k = (idx_vertex_lowest + j) % 3;
             int   l = (k + 1) % 3;
-            edge_id = edge_neighbor[k];
+            edge_id = edge_neighbor(k);
             a_id    = indices[k];
             a       = vertices + k;
             b_id    = indices[l];
@@ -147,7 +154,7 @@ static FacetSliceType slice_facet(
             FacetSliceType    result = FacetSliceType::Slicing;
             if (horizontal) {
                 // All three vertices are aligned with slice_z.
-                line_out->edge_type = IntersectionLine::FacetEdgeType::Horizontal;
+                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) {
@@ -165,19 +172,19 @@ static FacetSliceType slice_facet(
                 // in respect to the cutting plane).
                 result = third_below ? FacetSliceType::Slicing : FacetSliceType::Cutting;
                 if (third_below) {
-                    line_out->edge_type = IntersectionLine::FacetEdgeType::Top;
+                    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.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);
+            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;
         }
 
@@ -235,31 +242,31 @@ static FacetSliceType slice_facet(
     // 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;
+        line_out.edge_type  = IntersectionLine::FacetEdgeType::General;
+        line_out.a          = static_cast<const Point&>(points[1]);
+        line_out.b          = static_cast<const 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);
+        //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);
+        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) {
+        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)
+            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)
+            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 ?
+            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;
         }
@@ -269,40 +276,64 @@ static FacetSliceType slice_facet(
     return FacetSliceType::NoSlice;
 }
 
-static void slice_facet_at_zs(
+template<typename TransformVertex>
+void slice_facet_at_zs(
+    // Scaled or unscaled vertices. transform_vertex_fn may scale zs.
+    const std::vector<Vec3f>                         &mesh_vertices,
+    const TransformVertex                            &transform_vertex_fn,
     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)
+    const Vec3i                                      &facet_neighbors,
+    // Scaled or unscaled zs. If vertices have their zs scaled or transform_vertex_fn scales them, then zs have to be scaled as well.
+    const std::vector<float>                         &zs,
+    std::vector<IntersectionLines>                   &lines,
+    boost::mutex                                     &lines_mutex)
 {
-    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];
+    stl_vertex vertices[3] { transform_vertex_fn(mesh_vertices[indices(0)]), transform_vertex_fn(mesh_vertices[indices(1)]), transform_vertex_fn(mesh_vertices[indices(2)]) };
 
     // 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
+    auto min_layer = std::lower_bound(zs.begin(), zs.end(), min_z); // first layer whose slice_z is >= min_z
+    auto max_layer = std::upper_bound(min_layer, 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 &&
+        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);
+            boost::lock_guard<boost::mutex> l(lines_mutex);
+            lines[it - zs.begin()].emplace_back(il);
         }
     }
 }
 
+template<typename TransformVertex, typename ThrowOnCancel>
+inline std::vector<IntersectionLines> slice_make_lines(
+    const std::vector<stl_vertex>                   &vertices,
+    const TransformVertex                           &transform_vertex_fn,
+    const std::vector<stl_triangle_vertex_indices>  &indices,
+    const std::vector<Vec3i>                        &face_neighbors,
+    const std::vector<float>                        &zs,
+    const ThrowOnCancel                              throw_on_cancel_fn)
+{
+    std::vector<IntersectionLines>  lines(zs.size(), IntersectionLines());
+    boost::mutex                    lines_mutex;
+    tbb::parallel_for(
+        tbb::blocked_range<int>(0, int(indices.size())),
+        [&vertices, &transform_vertex_fn, &indices, &face_neighbors, &zs, &lines, &lines_mutex, throw_on_cancel_fn](const tbb::blocked_range<int> &range) {
+            for (int face_idx = range.begin(); face_idx < range.end(); ++ face_idx) {
+                if ((face_idx & 0x0ffff) == 0)
+                    throw_on_cancel_fn();
+                slice_facet_at_zs(vertices, transform_vertex_fn, indices[face_idx], face_neighbors[face_idx], zs, lines, lines_mutex);
+            }
+        }
+    );
+    return lines;
+}
+
 #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.
@@ -384,9 +415,9 @@ struct OpenPolyline {
     bool                    consumed;
 };
 
-// called by TriangleMeshSlicer::make_loops() to connect sliced triangles into closed loops and open polylines by the triangle connectivity.
+// called by 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)
+static void chain_lines_by_triangle_connectivity(IntersectionLines &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;
@@ -501,7 +532,7 @@ std::vector<OpenPolyline*> open_polylines_sorted(std::vector<OpenPolyline> &open
     return out;
 }
 
-// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices.
+// called by 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)
 {
@@ -599,7 +630,7 @@ static void chain_open_polylines_exact(std::vector<OpenPolyline> &open_polylines
     }
 }
 
-// called by TriangleMeshSlicer::make_loops() to connect remaining open polylines across shared triangle edges and vertices, 
+// called by 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)
@@ -707,8 +738,11 @@ static void chain_open_polylines_close_gaps(std::vector<OpenPolyline> &open_poly
     }
 }
 
-static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
+static Polygons make_loops(
+    // Lines will have their flags modified.
+    IntersectionLines   &lines)
 {
+    Polygons loops;
 #if 0
 //FIXME slice_facet() may create zero length edges due to rounding of doubles into coord_t.
 //#ifdef _DEBUG
@@ -736,7 +770,7 @@ static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
 #endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
 
     std::vector<OpenPolyline> open_polylines;
-    chain_lines_by_triangle_connectivity(lines, *loops, open_polylines);
+    chain_lines_by_triangle_connectivity(lines, loops, open_polylines);
 
 #ifdef SLIC3R_DEBUG_SLICE_PROCESSING
         {
@@ -752,8 +786,8 @@ static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
     // 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);
+    chain_open_polylines_exact(open_polylines, loops, false);
+    chain_open_polylines_exact(open_polylines, loops, true);
 
 #ifdef SLIC3R_DEBUG_SLICE_PROCESSING
     {
@@ -781,8 +815,8 @@ static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
     }
 #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);
+    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
@@ -800,6 +834,60 @@ static void make_loops(std::vector<IntersectionLine> &lines, Polygons* loops)
         svg.Close();
     }
 #endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
+
+    return loops;
+}
+
+template<typename ThrowOnCancel>
+static std::vector<Polygons> make_loops(
+    // Lines will have their flags modified.
+    std::vector<IntersectionLines> &lines, 
+    const MeshSlicingParams        &params, 
+    ThrowOnCancel                   throw_on_cancel)
+{
+    std::vector<Polygons> layers;
+    layers.resize(lines.size());
+    tbb::parallel_for(
+        tbb::blocked_range<size_t>(0, lines.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];
+                polygons = make_loops(lines[line_idx]);
+
+                auto this_mode = line_idx < params.slicing_mode_normal_below_layer ? params.mode_below : params.mode;
+                if (! polygons.empty()) {
+                    if (this_mode == MeshSlicingParams::SlicingMode::Positive) {
+                        // Reorient all loops to be CCW.
+                        for (Polygon& p : polygons)
+                            p.make_counter_clockwise();
+                    }
+                    else if (this_mode == MeshSlicingParams::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));
+                    }
+                }
+            }
+        }
+    );
+
+    return layers;
 }
 
 // Used to cut the mesh into two halves.
@@ -808,15 +896,11 @@ 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));
-    }
+    for (Polygon &loop : make_loops(lines))
+        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());
@@ -961,45 +1045,13 @@ static void make_expolygons(const Polygons &loops, const float closing_radius, c
         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,
+std::vector<Polygons> slice_mesh(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
     const MeshSlicingParams          &params,
-    std::vector<Polygons>           *layers,
-    throw_on_cancel_callback_type    throw_on_cancel) const
+    std::function<void()>             throw_on_cancel)
 {
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::slice";
+    BOOST_LOG_TRIVIAL(debug) << "slice_mesh to polygons";
 
     /*
        This method gets called with a list of unscaled Z coordinates and outputs
@@ -1027,72 +1079,28 @@ void TriangleMeshSlicer::slice(
         NOTE: this method accepts a vector of floats because the mesh coordinate
         type is float.
     */
+
+    std::vector<IntersectionLines> lines;
     
-    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)
+        std::vector<float> scaled_zs(zs);
+        for (float &z : scaled_zs)
             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);
-                }
-            }
-        );
+
+        std::vector<stl_vertex> v_scaled_shared(mesh.vertices);
+        for (stl_vertex &v : v_scaled_shared)
+            v *= float(1. / SCALING_FACTOR);
+
+        std::vector<Vec3i> facets_edges = create_face_neighbors_index(mesh);
+        lines = params.trafo.matrix() == Transform3f::Identity().matrix() ?
+            slice_make_lines(v_scaled_shared, [](const Vec3f &p) { return p; },  mesh.indices, facets_edges, scaled_zs, throw_on_cancel) :
+            slice_make_lines(v_scaled_shared, [&params](const Vec3f &p) { return params.trafo * p; }, mesh.indices, facets_edges, scaled_zs, throw_on_cancel);
+        throw_on_cancel();
     }
-    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";
+    std::vector<Polygons> layers = make_loops(lines, params, throw_on_cancel);
 
 #ifdef SLIC3R_DEBUG
     {
@@ -1126,65 +1134,104 @@ void TriangleMeshSlicer::slice(
         ++ iRun;
     }
 #endif
+
+    return layers;
 }
 
-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<ExPolygons> slice_mesh_ex(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
+    const MeshSlicingParamsEx        &params,
+    std::function<void()>             throw_on_cancel)
 {
     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);
+        if (params.mode == MeshSlicingParams::SlicingMode::PositiveLargestContour)
+            slicing_params.mode = MeshSlicingParams::SlicingMode::Positive;
+        if (params.mode_below == MeshSlicingParams::SlicingMode::PositiveLargestContour)
+            slicing_params.mode_below = MeshSlicingParams::SlicingMode::Positive;
+        layers_p = slice_mesh(mesh, zs, slicing_params, throw_on_cancel);
     }
     
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - start";
-    layers->resize(z.size());
+    BOOST_LOG_TRIVIAL(debug) << "slice_mesh make_expolygons in parallel - start";
+    std::vector<ExPolygons> layers(layers_p.size(), ExPolygons{});
     tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, z.size()),
-        [&layers_p, &params, layers, throw_on_cancel]
+        tbb::blocked_range<size_t>(0, layers_p.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];
+                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)
+                if (this_mode == MeshSlicingParams::SlicingMode::PositiveLargestContour)
                     keep_largest_contour_only(expolygons);
             }
         });
-    BOOST_LOG_TRIVIAL(debug) << "TriangleMeshSlicer::make_expolygons in parallel - end";
+    BOOST_LOG_TRIVIAL(debug) << "slice_mesh make_expolygons in parallel - end";
+
+    return layers;
 }
 
-static void triangulate_slice(indexed_triangle_set &its, IntersectionLines &lines, const std::vector<int> &slice_vertices, float z)
+static void triangulate_slice(
+    indexed_triangle_set    &its, 
+    IntersectionLines       &lines, 
+    std::vector<int>        &slice_vertices,
+    // Vertices of the original (unsliced) mesh. Newly added vertices are those on the slice.
+    int                      num_original_vertices,
+    // Z height of the slice.
+    float                    z)
 {
-//    BOOST_LOG_TRIVIAL(trace) << "TriangleMeshSlicer::cut - triangulating the cut";
+    sort_remove_duplicates(slice_vertices);
 
-    ExPolygons section   = make_expolygons_simple(lines);
-    Pointf3s   triangles = triangulate_expolygons_3d(section, z, true);
+    // 1) Create map of the slice vertices from positions to mesh indices.
+    // As the caller will likely add duplicate points when intersecting triangle edges, there will be duplicates.
     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()); });
+            return l.first.x() < r.first.x() || 
+                   (l.first.x() == r.first.x() && (l.first.y() < r.first.y() || 
+                                                   (l.first.y() == r.first.y() && l.second < r.second))); });
+
+    // 2) Discover duplicate points on the slice. Remap duplicate vertices to a vertex with a lowest index.
+    {
+        std::vector<int> map_duplicate_vertex(int(its.vertices.size()) - num_original_vertices, -1);
+        int i = 0;
+        for (; i < int(map_vertex_to_index.size()); ++ i) {
+            const Vec2f &ipos = map_vertex_to_index[i].first;
+            const int    iidx = map_vertex_to_index[i].second;
+            if (iidx >= num_original_vertices)
+                // map to itself
+                map_duplicate_vertex[iidx - num_original_vertices] = iidx;
+            int j = i;
+            for (++ j; j < int(map_vertex_to_index.size()) && ipos.x() == map_vertex_to_index[j].first.x() && ipos.y() == map_vertex_to_index[j].first.y(); ++ j) {
+                const int jidx = map_vertex_to_index[j].second;
+                assert(jidx >= num_original_vertices);
+                if (jidx >= num_original_vertices)
+                    // map to the first vertex
+                    map_duplicate_vertex[jidx - num_original_vertices] = iidx;
+            }
+        }
+        for (stl_triangle_vertex_indices &f : its.indices)
+            for (i = 0; i < 3; ++ i)
+                if (f(i) >= num_original_vertices)
+                    f(i) = map_duplicate_vertex[f(i) - num_original_vertices];
+    }
+
     size_t idx_vertex_new_first = its.vertices.size();
+    Pointf3s triangles = triangulate_expolygons_3d(make_expolygons_simple(lines), z, true);
     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>();
+            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;
@@ -1204,48 +1251,64 @@ static void triangulate_slice(indexed_triangle_set &its, IntersectionLines &line
             }
             facet(j) = idx;
         }
-        its.indices.emplace_back(facet);
+        if (facet(0) != facet(1) && facet(0) != facet(2) && facet(1) != facet(2))
+            its.indices.emplace_back(facet);
     }
 
+    // Remove vertices, which are not referenced by any face.
     its_compactify_vertices(its);
-    its_remove_degenerate_faces(its);
+
+    // Degenerate faces should not be created.
+    // its_remove_degenerate_faces(its);
 }
 
-void TriangleMeshSlicer::cut(float z, indexed_triangle_set *upper, indexed_triangle_set *lower) const
+void cut_mesh(const indexed_triangle_set &mesh, float z, indexed_triangle_set *upper, indexed_triangle_set *lower)
 {
     assert(upper || lower);
     if (upper == nullptr && lower == nullptr)
         return;
 
+    BOOST_LOG_TRIVIAL(trace) << "cut_mesh - slicing object";
+
     if (upper) {
         upper->clear();
-        upper->vertices = m_its->vertices;
-        upper->indices.reserve(m_its->indices.size());
+        upper->vertices = mesh.vertices;
+        upper->indices.reserve(mesh.indices.size());
     }
+
     if (lower) {
         lower->clear();
-        lower->vertices = m_its->vertices;
-        lower->indices.reserve(m_its->indices.size());
+        lower->vertices = mesh.vertices;
+        lower->indices.reserve(mesh.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;
+    IntersectionLines  upper_lines, lower_lines;
+    std::vector<int>   upper_slice_vertices, lower_slice_vertices;
+    std::vector<Vec3i> facets_edges = create_face_neighbors_index(mesh);
 
-    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]] };
+    for (int facet_idx = 0; facet_idx < int(mesh.indices.size()); ++ facet_idx) {
+        const stl_triangle_vertex_indices &facet = mesh.indices[facet_idx];
+        Vec3f vertices[3] { mesh.vertices[facet(0)], mesh.vertices[facet(1)], mesh.vertices[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);
+        FacetSliceType   slice_type = FacetSliceType::NoSlice;
+        if (z > min_z - EPSILON && z < max_z + EPSILON) {
+            Vec3f vertices_scaled[3];
+            for (int i = 0; i < 3; ++ i) {
+                const Vec3f &src = vertices[i];
+                Vec3f       &dst = vertices_scaled[i];
+                dst.x() = scale_(src.x());
+                dst.y() = scale_(src.y());
+                dst.z() = src.z();
+            }
+            slice_type = slice_facet(z, vertices_scaled, mesh.indices[facet_idx], facets_edges[facet_idx], 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) {
@@ -1274,6 +1337,8 @@ void TriangleMeshSlicer::cut(float z, indexed_triangle_set *upper, indexed_trian
             // Facet is cut by the slicing plane.
             assert(slice_type == FacetSliceType::Slicing);
             assert(line.edge_type == IntersectionLine::FacetEdgeType::General);
+            assert(line.edge_a_id != -1);
+            assert(line.edge_b_id != -1);
 
             // look for the vertex on whose side of the slicing plane there are no other vertices
             int isolated_vertex = 
@@ -1282,6 +1347,15 @@ void TriangleMeshSlicer::cut(float z, indexed_triangle_set *upper, indexed_trian
             
             // get vertices starting from the isolated one
             int iv = isolated_vertex;
+            stl_vertex v0v1, v2v0;
+            assert(facets_edges[facet_idx](iv) == line.edge_a_id ||facets_edges[facet_idx](iv) == line.edge_b_id);
+            if (facets_edges[facet_idx](iv) == line.edge_a_id) {
+                v0v1 = to_3d(unscaled<float>(line.a), z);
+                v2v0 = to_3d(unscaled<float>(line.b), z);
+            } else {
+                v0v1 = to_3d(unscaled<float>(line.b), z);
+                v2v0 = to_3d(unscaled<float>(line.a), z);
+            }
             const stl_vertex &v0  = vertices[iv];
             const int         iv0 = facet[iv];
             if (++ iv == 3)
@@ -1294,70 +1368,51 @@ void TriangleMeshSlicer::cut(float z, indexed_triangle_set *upper, indexed_trian
             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) {
+            auto new_vertex = [upper, lower, &upper_slice_vertices, &lower_slice_vertices](const Vec3f &a, const int ia, const Vec3f &b, const int ib, const Vec3f &c) {
                 int iupper, ilower;
-                if (t <= 0.f)
+                if (c == a)
                     iupper = ilower = ia;
-                else if (t >= 1.f)
+                else if (c == b)
                     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);
-                        }
+                    // 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);
-                }
+            auto [iv0v1_upper, iv0v1_lower] = new_vertex(v1, iv1, v0, iv0, v0v1);
+            auto [iv2v0_upper, iv2v0_lower] = new_vertex(v2, iv2, v0, iv0, v2v0);
+            auto new_face = [](indexed_triangle_set *its, int i, int j, int k) {
+                if (its != nullptr && i != j && i != k && j != k)
+                    its->indices.emplace_back(i, j, k);
+            };
+
+            if (v0.z() > z) {
+                new_face(upper, iv0, iv0v1_upper, iv2v0_upper);
+                new_face(lower, iv1, iv2, iv0v1_lower);
+                new_face(lower, 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);
+                new_face(upper, iv1, iv2, iv0v1_upper);
+                new_face(upper, iv2, iv2v0_upper, iv0v1_upper);
+                new_face(lower, iv0, iv0v1_lower, iv2v0_lower);
             }
         }
     }
     
     if (upper != nullptr)
-        triangulate_slice(*upper, upper_lines, upper_slice_vertices, z);
+        triangulate_slice(*upper, upper_lines, upper_slice_vertices, int(mesh.vertices.size()), 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);
+        triangulate_slice(*lower, lower_lines, lower_slice_vertices, int(mesh.vertices.size()), z);
 }
 
 }
diff --git a/src/libslic3r/TriangleMeshSlicer.hpp b/src/libslic3r/TriangleMeshSlicer.hpp
index 258ffa04c..b4914f88b 100644
--- a/src/libslic3r/TriangleMeshSlicer.hpp
+++ b/src/libslic3r/TriangleMeshSlicer.hpp
@@ -1,42 +1,36 @@
 #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
 {
+    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,
+    };
+
     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 };
+    // Transforming faces during the slicing.
+    Transform3f   trafo { Transform3f::Identity() };
 };
 
-struct MeshSlicingParamsExtended : public MeshSlicingParams
+struct MeshSlicingParamsEx : public MeshSlicingParams
 {
     // Morphological closing operation when creating output expolygons.
     float         closing_radius { 0 };
@@ -45,96 +39,44 @@ struct MeshSlicingParamsExtended : public MeshSlicingParams
     // Resolution for contour simplification.
     // 0 = don't simplify.
     double        resolution { 0 };
-    // Transformation of the object owning the ModelVolume.
-//    Transform3d         object_trafo;
 };
 
-class TriangleMeshSlicer
+std::vector<Polygons>           slice_mesh(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
+    const MeshSlicingParams          &params,
+    std::function<void()>             throw_on_cancel = []{});
+
+std::vector<ExPolygons>         slice_mesh_ex(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
+    const MeshSlicingParamsEx        &params,
+    std::function<void()>             throw_on_cancel = []{});
+
+inline std::vector<ExPolygons>  slice_mesh_ex(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
+    std::function<void()>             throw_on_cancel = []{})
 {
-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);
-    }
+    return slice_mesh_ex(mesh, zs, MeshSlicingParamsEx{}, throw_on_cancel);
 }
 
-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 = []{})
+inline std::vector<ExPolygons>  slice_mesh_ex(
+    const indexed_triangle_set       &mesh,
+    const std::vector<float>         &zs,
+    float                             closing_radius,
+    std::function<void()>             throw_on_cancel = []{})
 {
-    MeshSlicingParamsExtended params;
+    MeshSlicingParamsEx params;
     params.closing_radius = closing_radius;
-    slice_mesh(mesh, z, params, layers);
+    return slice_mesh_ex(mesh, zs, params, throw_on_cancel);
 }
 
-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);
-}
+void                            cut_mesh(
+    const indexed_triangle_set      &mesh,
+    float                            z,
+    indexed_triangle_set            *upper,
+    indexed_triangle_set            *lower);
 
 }