diff --git a/sandboxes/CMakeLists.txt b/sandboxes/CMakeLists.txt
index 23c15f0899..81ea94a291 100644
--- a/sandboxes/CMakeLists.txt
+++ b/sandboxes/CMakeLists.txt
@@ -1,5 +1,6 @@
 #add_subdirectory(slasupporttree)
 #add_subdirectory(openvdb)
-add_subdirectory(meshboolean)
-add_subdirectory(opencsg)
+# add_subdirectory(meshboolean)
+add_subdirectory(its_neighbor_index)
+# add_subdirectory(opencsg)
 #add_subdirectory(aabb-evaluation)
\ No newline at end of file
diff --git a/sandboxes/its_neighbor_index/CMakeLists.txt b/sandboxes/its_neighbor_index/CMakeLists.txt
new file mode 100644
index 0000000000..e289399482
--- /dev/null
+++ b/sandboxes/its_neighbor_index/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_executable(its_neighbor_index main.cpp ItsNeighborIndex.cpp ItsNeighborIndex.hpp)
+
+target_link_libraries(its_neighbor_index libslic3r admesh)
+
+if (WIN32)
+    prusaslicer_copy_dlls(its_neighbor_index)
+endif()
diff --git a/sandboxes/its_neighbor_index/ItsNeighborIndex.cpp b/sandboxes/its_neighbor_index/ItsNeighborIndex.cpp
new file mode 100644
index 0000000000..92c0f57c33
--- /dev/null
+++ b/sandboxes/its_neighbor_index/ItsNeighborIndex.cpp
@@ -0,0 +1,613 @@
+#include <iostream>
+#include <vector>
+#include <unordered_map>
+#include <map>
+
+#include "ItsNeighborIndex.hpp"
+#include "libslic3r/Execution/ExecutionTBB.hpp"
+#include "libslic3r/Execution/ExecutionSeq.hpp"
+
+#include "tbb/parallel_sort.h"
+
+namespace Slic3r {
+
+FaceNeighborIndex its_create_neighbors_index_1(const indexed_triangle_set &its)
+{
+    // Just to be clear what type of object are we referencing
+    using FaceID = size_t;
+    using VertexID = uint64_t;
+    using EdgeID = uint64_t;
+
+    constexpr auto UNASSIGNED = std::numeric_limits<FaceID>::max();
+
+    struct Edge // Will contain IDs of the two facets touching this edge
+    {
+        FaceID first, second;
+        Edge() : first{UNASSIGNED}, second{UNASSIGNED} {}
+        void   assign(FaceID fid)
+        {
+            first == UNASSIGNED ? first = fid : second = fid;
+        }
+    };
+
+    // All vertex IDs will fit into this number of bits. (Used for hashing)
+    const int max_vertex_id_bits = std::ceil(std::log2(its.vertices.size()));
+    assert(max_vertex_id_bits <= 32);
+
+    std::unordered_map< EdgeID, Edge> edge_index;
+
+    // Edge id is constructed by concatenating two vertex ids, starting with
+    // the lowest in MSB
+    auto hash = [max_vertex_id_bits] (VertexID a, VertexID b) {
+        if (a > b) std::swap(a, b);
+        return (a << max_vertex_id_bits) + b;
+    };
+
+    // Go through all edges of all facets and mark the facets touching each edge
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        edge_index[e1].assign(face_id);
+        edge_index[e2].assign(face_id);
+        edge_index[e3].assign(face_id);
+    }
+
+    FaceNeighborIndex index(its.indices.size());
+
+    // Now collect the neighbors for each facet into the final index
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        const Edge &neighs1 = edge_index[e1];
+        const Edge &neighs2 = edge_index[e2];
+        const Edge &neighs3 = edge_index[e3];
+
+        std::array<size_t, 3> &neighs = index[face_id];
+        neighs[0] = neighs1.first == face_id ? neighs1.second : neighs1.first;
+        neighs[1] = neighs2.first == face_id ? neighs2.second : neighs2.first;
+        neighs[2] = neighs3.first == face_id ? neighs3.second : neighs3.first;
+    }
+
+    return index;
+}
+
+std::vector<Vec3i> its_create_neighbors_index_2(const indexed_triangle_set &its)
+{
+    std::vector<Vec3i> out(its.indices.size(), Vec3i(-1, -1, -1));
+
+    // Create a mapping from triangle edge into face.
+    struct EdgeToFace {
+        // Index of the 1st vertex of the triangle edge. vertex_low <= vertex_high.
+        int  vertex_low;
+        // Index of the 2nd vertex of the triangle edge.
+        int  vertex_high;
+        // Index of a triangular face.
+        int  face;
+        // Index of edge in the face, starting with 1. Negative indices if the edge was stored reverse in (vertex_low, vertex_high).
+        int  face_edge;
+        bool operator==(const EdgeToFace &other) const { return vertex_low == other.vertex_low && vertex_high == other.vertex_high; }
+        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(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  = 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;
+            if (e2f.vertex_low > e2f.vertex_high) {
+                // Sort the vertices
+                std::swap(e2f.vertex_low, e2f.vertex_high);
+                // and make the face_edge negative to indicate a flipped edge.
+                e2f.face_edge = - e2f.face_edge;
+            }
+        }
+
+    std::sort(edges_map.begin(), edges_map.end());
+
+    // Assign a unique common edge id to touching triangle edges.
+    int num_edges = 0;
+    for (size_t i = 0; i < edges_map.size(); ++ i) {
+        EdgeToFace &edge_i = edges_map[i];
+        if (edge_i.face == -1)
+            // This edge has been connected to some neighbor already.
+            continue;
+        // Unconnected edge. Find its neighbor with the correct orientation.
+        size_t j;
+        bool found = false;
+        for (j = i + 1; j < edges_map.size() && edge_i == edges_map[j]; ++ j)
+            if (edge_i.face_edge * edges_map[j].face_edge < 0 && edges_map[j].face != -1) {
+                // Faces touching with opposite oriented edges and none of the edges is connected yet.
+                found = true;
+                break;
+            }
+        if (! found) {
+            //FIXME Vojtech: Trying to find an edge with equal orientation. This smells.
+            // admesh can assign the same edge ID to more than two facets (which is
+            // still topologically correct), so we have to search for a duplicate of
+            // this edge too in case it was already seen in this orientation
+            for (j = i + 1; j < edges_map.size() && edge_i == edges_map[j]; ++ j)
+                if (edges_map[j].face != -1) {
+                    // Faces touching with equally oriented edges and none of the edges is connected yet.
+                    found = true;
+                    break;
+                }
+        }
+        // Assign an edge index to the 1st face.
+        //        out[edge_i.face](std::abs(edge_i.face_edge) - 1) = num_edges;
+        if (found) {
+            EdgeToFace &edge_j = edges_map[j];
+            out[edge_i.face](std::abs(edge_i.face_edge) - 1) = edge_j.face;
+            out[edge_j.face](std::abs(edge_j.face_edge) - 1) = edge_i.face;
+            // Mark the edge as connected.
+            edge_j.face = -1;
+        }
+        ++ num_edges;
+    }
+
+    return out;
+}
+
+std::vector<Vec3i> its_create_neighbors_index_3(const indexed_triangle_set &its)
+{
+    std::vector<Vec3i> out(its.indices.size(), Vec3i(-1, -1, -1));
+
+    // Create a mapping from triangle edge into face.
+    struct EdgeToFace {
+        // Index of the 1st vertex of the triangle edge. vertex_low <= vertex_high.
+        int  vertex_low;
+        // Index of the 2nd vertex of the triangle edge.
+        int  vertex_high;
+        // Index of a triangular face.
+        int  face;
+        // Index of edge in the face, starting with 1. Negative indices if the edge was stored reverse in (vertex_low, vertex_high).
+        int  face_edge;
+        bool operator==(const EdgeToFace &other) const { return vertex_low == other.vertex_low && vertex_high == other.vertex_high; }
+        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(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  = 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;
+            if (e2f.vertex_low > e2f.vertex_high) {
+                // Sort the vertices
+                std::swap(e2f.vertex_low, e2f.vertex_high);
+                // and make the face_edge negative to indicate a flipped edge.
+                e2f.face_edge = - e2f.face_edge;
+            }
+        }
+
+    tbb::parallel_sort(edges_map.begin(), edges_map.end());
+
+    // Assign a unique common edge id to touching triangle edges.
+    int num_edges = 0;
+    for (size_t i = 0; i < edges_map.size(); ++ i) {
+        EdgeToFace &edge_i = edges_map[i];
+        if (edge_i.face == -1)
+            // This edge has been connected to some neighbor already.
+            continue;
+        // Unconnected edge. Find its neighbor with the correct orientation.
+        size_t j;
+        bool found = false;
+        for (j = i + 1; j < edges_map.size() && edge_i == edges_map[j]; ++ j)
+            if (edge_i.face_edge * edges_map[j].face_edge < 0 && edges_map[j].face != -1) {
+                // Faces touching with opposite oriented edges and none of the edges is connected yet.
+                found = true;
+                break;
+            }
+        if (! found) {
+            //FIXME Vojtech: Trying to find an edge with equal orientation. This smells.
+            // admesh can assign the same edge ID to more than two facets (which is
+            // still topologically correct), so we have to search for a duplicate of
+            // this edge too in case it was already seen in this orientation
+            for (j = i + 1; j < edges_map.size() && edge_i == edges_map[j]; ++ j)
+                if (edges_map[j].face != -1) {
+                    // Faces touching with equally oriented edges and none of the edges is connected yet.
+                    found = true;
+                    break;
+                }
+        }
+        // Assign an edge index to the 1st face.
+        //        out[edge_i.face](std::abs(edge_i.face_edge) - 1) = num_edges;
+        if (found) {
+            EdgeToFace &edge_j = edges_map[j];
+            out[edge_i.face](std::abs(edge_i.face_edge) - 1) = edge_j.face;
+            out[edge_j.face](std::abs(edge_j.face_edge) - 1) = edge_i.face;
+            // Mark the edge as connected.
+            edge_j.face = -1;
+        }
+        ++ num_edges;
+    }
+
+    return out;
+}
+
+FaceNeighborIndex its_create_neighbors_index_4(const indexed_triangle_set &its)
+{
+    // Just to be clear what type of object are we referencing
+    using FaceID = size_t;
+    using VertexID = uint64_t;
+    using EdgeID = uint64_t;
+
+    constexpr auto UNASSIGNED = std::numeric_limits<FaceID>::max();
+
+    struct Edge // Will contain IDs of the two facets touching this edge
+    {
+        FaceID first, second;
+        Edge() : first{UNASSIGNED}, second{UNASSIGNED} {}
+        void   assign(FaceID fid)
+        {
+            first == UNASSIGNED ? first = fid : second = fid;
+        }
+    };
+
+    Benchmark bm;
+    bm.start();
+
+    // All vertex IDs will fit into this number of bits. (Used for hashing)
+    //    const int max_vertex_id_bits = std::ceil(std::log2(its.vertices.size()));
+    //    assert(max_vertex_id_bits <= 32);
+
+    const uint64_t Vn  = its.vertices.size();
+    //    const uint64_t Fn  = 3 * its.indices.size();
+    //    double MaxQ = double(Vn) * (Vn + 1) / Fn;
+    //    const uint64_t Nq = MaxQ < 0 ? 0 : std::ceil(std::log2(MaxQ));
+    //    const uint64_t Nr = std::ceil(std::log2(std::min(Vn * (Vn + 1), Fn)));
+    //    const uint64_t Nfn = std::ceil(std::log2(Fn));
+
+    ////    const uint64_t max_edge_ids = (uint64_t(1) << (Nq + Nr));
+    //    const uint64_t max_edge_ids = MaxQ * Fn + (std::min(Vn * (Vn + 1), Fn)); //(uint64_t(1) << Nfn);
+    const uint64_t Fn  = 3 * its.indices.size();
+    std::vector< Edge > edge_index(3 * Fn);
+
+    // Edge id is constructed by concatenating two vertex ids, starting with
+    // the lowest in MSB
+    auto hash = [Vn, Fn /*, Nr*/] (VertexID a, VertexID b) {
+        if (a > b) std::swap(a, b);
+
+        uint64_t C = Vn * a + b;
+        uint64_t Q = C / Fn, R = C % Fn;
+
+        return Q * Fn + R;
+    };
+
+    // Go through all edges of all facets and mark the facets touching each edge
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        edge_index[e1].assign(face_id);
+        edge_index[e2].assign(face_id);
+        edge_index[e3].assign(face_id);
+    }
+
+    FaceNeighborIndex index(its.indices.size());
+
+    // Now collect the neighbors for each facet into the final index
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        const Edge &neighs1 = edge_index[e1];
+        const Edge &neighs2 = edge_index[e2];
+        const Edge &neighs3 = edge_index[e3];
+
+        std::array<size_t, 3> &neighs = index[face_id];
+        neighs[0] = neighs1.first == face_id ? neighs1.second : neighs1.first;
+        neighs[1] = neighs2.first == face_id ? neighs2.second : neighs2.first;
+        neighs[2] = neighs3.first == face_id ? neighs3.second : neighs3.first;
+    }
+
+    bm.stop();
+
+    std::cout << "Creating neighbor index took: " << bm.getElapsedSec() << " seconds." << std::endl;
+
+    return index;
+}
+
+// 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_vertex_faces_index(const indexed_triangle_set &its)
+{
+    std::vector<std::vector<size_t>> index;
+
+    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;
+}
+
+//static int get_vertex_index(size_t vertex_index, const stl_triangle_vertex_indices &triangle_indices) {
+//    if (vertex_index == triangle_indices[0]) return 0;
+//    if (vertex_index == triangle_indices[1]) return 1;
+//    if (vertex_index == triangle_indices[2]) return 2;
+//    return -1;
+//}
+
+//static Vec2crd get_edge_indices(int edge_index, const stl_triangle_vertex_indices &triangle_indices)
+//{
+//    int next_edge_index = (edge_index == 2) ? 0 : edge_index + 1;
+//    coord_t vi0             = triangle_indices[edge_index];
+//    coord_t vi1             = triangle_indices[next_edge_index];
+//    return Vec2crd(vi0, vi1);
+//}
+
+static std::vector<std::vector<size_t>> create_vertex_faces_index(
+    const std::vector<stl_triangle_vertex_indices>& indices, size_t count_vertices)
+{
+    if (count_vertices == 0) return {};
+    std::vector<std::vector<size_t>> index;
+    size_t res = indices.size() / count_vertices;
+    index.assign(count_vertices, reserve_vector<size_t>(res));
+    for (size_t fi = 0; fi < indices.size(); ++fi) {
+        auto &face = indices[fi];
+        index[face(0)].emplace_back(fi);
+        index[face(1)].emplace_back(fi);
+        index[face(2)].emplace_back(fi);
+    }
+    return index;
+}
+
+std::vector<Vec3crd> its_create_neighbors_index_5(const indexed_triangle_set &its)
+{
+    const std::vector<stl_triangle_vertex_indices> &indices = its.indices;
+    size_t vertices_size = its.vertices.size();
+
+    if (indices.empty() || vertices_size == 0) return {};
+    std::vector<std::vector<size_t>> vertex_triangles = create_vertex_faces_index(indices, vertices_size);
+    coord_t              no_value = -1;
+    std::vector<Vec3crd> neighbors(indices.size(), Vec3crd(no_value, no_value, no_value));
+    for (const stl_triangle_vertex_indices& triangle_indices : indices) {
+        coord_t index = &triangle_indices - &indices.front();
+        Vec3crd& neighbor = neighbors[index];
+        for (int edge_index = 0; edge_index < 3; ++edge_index) {
+            // check if done
+            coord_t& neighbor_edge = neighbor[edge_index];
+            if (neighbor_edge != no_value) continue;
+            Vec2crd edge_indices = get_edge_indices(edge_index, triangle_indices);
+            // IMPROVE: use same vector for 2 sides of triangle
+            const std::vector<size_t> &faces = vertex_triangles[edge_indices[0]];
+            for (const size_t &face : faces) {
+                if (face <= index) continue;
+                const stl_triangle_vertex_indices &face_indices = indices[face];
+                int vertex_index = get_vertex_index(edge_indices[1], face_indices);
+                // NOT Contain second vertex?
+                if (vertex_index < 0) continue;
+                // Has NOT oposit direction?
+                if (edge_indices[0] != face_indices[(vertex_index + 1) % 3]) continue;
+                neighbor_edge = face;
+                neighbors[face][vertex_index] = index;
+                break;
+            }
+            // must be paired
+            assert(neighbor_edge != no_value);
+        }
+    }
+
+    return neighbors;
+}
+
+std::vector<std::array<size_t, 3>> its_create_neighbors_index_6(const indexed_triangle_set &its)
+{
+    constexpr auto UNASSIGNED_EDGE = std::numeric_limits<uint64_t>::max();
+    constexpr auto UNASSIGNED_FACE = std::numeric_limits<size_t>::max();
+    struct Edge
+    {
+        uint64_t id      = UNASSIGNED_EDGE;
+        size_t   face_id = UNASSIGNED_FACE;
+        bool operator < (const Edge &e) const { return id < e.id; }
+    };
+
+    const size_t facenum = its.indices.size();
+
+    // All vertex IDs will fit into this number of bits. (Used for hashing)
+    const int max_vertex_id_bits = std::ceil(std::log2(its.vertices.size()));
+    assert(max_vertex_id_bits <= 32);
+
+    // Edge id is constructed by concatenating two vertex ids, starting with
+    // the lowest in MSB
+    auto hash = [max_vertex_id_bits] (uint64_t a, uint64_t b) {
+        if (a > b) std::swap(a, b);
+        return (a << max_vertex_id_bits) + b;
+    };
+
+    std::vector<Edge> edge_map(3 * facenum);
+
+    // Go through all edges of all facets and mark the facets touching each edge
+    for (size_t face_id = 0; face_id < facenum; ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        edge_map[face_id * 3] = {hash(face(0), face(1)), face_id};
+        edge_map[face_id * 3 + 1] = {hash(face(1), face(2)), face_id};
+        edge_map[face_id * 3 + 2] = {hash(face(2), face(0)), face_id};
+    }
+
+    std::sort(edge_map.begin(), edge_map.end());
+
+    std::vector<std::array<size_t, 3>> out(facenum, {UNASSIGNED_FACE, UNASSIGNED_FACE, UNASSIGNED_FACE});
+
+    auto add_neighbor = [](std::array<size_t, 3> &slot, size_t face_id) {
+        if (slot[0] == UNASSIGNED_FACE) { slot[0] = face_id; return; }
+        if (slot[1] == UNASSIGNED_FACE) { slot[1] = face_id; return; }
+        if (slot[2] == UNASSIGNED_FACE) { slot[2] = face_id; return; }
+    };
+
+    for (auto it = edge_map.begin(); it != edge_map.end();) {
+        size_t face_id = it->face_id;
+        uint64_t edge_id = it->id;
+
+        while (++it != edge_map.end() &&  (it->id == edge_id)) {
+            size_t other_face_id = it->face_id;
+            add_neighbor(out[other_face_id], face_id);
+            add_neighbor(out[face_id], other_face_id);
+        }
+    }
+
+    return out;
+}
+
+
+std::vector<std::array<size_t, 3>> its_create_neighbors_index_7(const indexed_triangle_set &its)
+{
+    constexpr auto UNASSIGNED_EDGE = std::numeric_limits<uint64_t>::max();
+    constexpr auto UNASSIGNED_FACE = std::numeric_limits<size_t>::max();
+    struct Edge
+    {
+        uint64_t id      = UNASSIGNED_EDGE;
+        size_t   face_id = UNASSIGNED_FACE;
+        bool operator < (const Edge &e) const { return id < e.id; }
+    };
+
+    const size_t facenum = its.indices.size();
+
+    // All vertex IDs will fit into this number of bits. (Used for hashing)
+    const int max_vertex_id_bits = std::ceil(std::log2(its.vertices.size()));
+    assert(max_vertex_id_bits <= 32);
+
+    // Edge id is constructed by concatenating two vertex ids, starting with
+    // the lowest in MSB
+    auto hash = [max_vertex_id_bits] (uint64_t a, uint64_t b) {
+        if (a > b) std::swap(a, b);
+        return (a << max_vertex_id_bits) + b;
+    };
+
+    std::vector<Edge> edge_map(3 * facenum);
+
+    // Go through all edges of all facets and mark the facets touching each edge
+    for (size_t face_id = 0; face_id < facenum; ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        edge_map[face_id * 3] = {hash(face(0), face(1)), face_id};
+        edge_map[face_id * 3 + 1] = {hash(face(1), face(2)), face_id};
+        edge_map[face_id * 3 + 2] = {hash(face(2), face(0)), face_id};
+    }
+
+    tbb::parallel_sort(edge_map.begin(), edge_map.end());
+
+    std::vector<std::array<size_t, 3>> out(facenum, {UNASSIGNED_FACE, UNASSIGNED_FACE, UNASSIGNED_FACE});
+
+    auto add_neighbor = [](std::array<size_t, 3> &slot, size_t face_id) {
+        if (slot[0] == UNASSIGNED_FACE) { slot[0] = face_id; return; }
+        if (slot[1] == UNASSIGNED_FACE) { slot[1] = face_id; return; }
+        if (slot[2] == UNASSIGNED_FACE) { slot[2] = face_id; return; }
+    };
+
+    for (auto it = edge_map.begin(); it != edge_map.end();) {
+        size_t face_id = it->face_id;
+        uint64_t edge_id = it->id;
+
+        while (++it != edge_map.end() &&  (it->id == edge_id)) {
+            size_t other_face_id = it->face_id;
+            add_neighbor(out[other_face_id], face_id);
+            add_neighbor(out[face_id], other_face_id);
+        }
+    }
+
+    return out;
+}
+
+FaceNeighborIndex its_create_neighbors_index_8(const indexed_triangle_set &its)
+{
+    // Just to be clear what type of object are we referencing
+    using FaceID = size_t;
+    using VertexID = uint64_t;
+    using EdgeID = uint64_t;
+
+    constexpr auto UNASSIGNED = std::numeric_limits<FaceID>::max();
+
+    struct Edge // Will contain IDs of the two facets touching this edge
+    {
+        FaceID first, second;
+        Edge() : first{UNASSIGNED}, second{UNASSIGNED} {}
+        void   assign(FaceID fid)
+        {
+            first == UNASSIGNED ? first = fid : second = fid;
+        }
+    };
+
+    // All vertex IDs will fit into this number of bits. (Used for hashing)
+    const int max_vertex_id_bits = std::ceil(std::log2(its.vertices.size()));
+    assert(max_vertex_id_bits <= 32);
+
+    std::map< EdgeID, Edge > edge_index;
+
+    // Edge id is constructed by concatenating two vertex ids, starting with
+    // the lowest in MSB
+    auto hash = [max_vertex_id_bits] (VertexID a, VertexID b) {
+        if (a > b) std::swap(a, b);
+        return (a << max_vertex_id_bits) + b;
+    };
+
+    // Go through all edges of all facets and mark the facets touching each edge
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        edge_index[e1].assign(face_id);
+        edge_index[e2].assign(face_id);
+        edge_index[e3].assign(face_id);
+    }
+
+    FaceNeighborIndex index(its.indices.size());
+
+    // Now collect the neighbors for each facet into the final index
+    for (size_t face_id = 0; face_id < its.indices.size(); ++face_id) {
+        const Vec3i &face = its.indices[face_id];
+
+        EdgeID e1 = hash(face(0), face(1)), e2 = hash(face(1), face(2)),
+               e3 = hash(face(2), face(0));
+
+        const Edge &neighs1 = edge_index[e1];
+        const Edge &neighs2 = edge_index[e2];
+        const Edge &neighs3 = edge_index[e3];
+
+        std::array<size_t, 3> &neighs = index[face_id];
+        neighs[0] = neighs1.first == face_id ? neighs1.second : neighs1.first;
+        neighs[1] = neighs2.first == face_id ? neighs2.second : neighs2.first;
+        neighs[2] = neighs3.first == face_id ? neighs3.second : neighs3.first;
+    }
+
+    return index;
+}
+
+std::vector<Vec3crd> its_create_neighbors_index_9(const indexed_triangle_set &its)
+{
+    return create_neighbors_index(ex_seq, its);
+}
+
+std::vector<Vec3i> its_create_neighbors_index_10(const indexed_triangle_set &its)
+{
+    return create_neighbors_index(ex_tbb, its);
+}
+
+} // namespace Slic3r
diff --git a/sandboxes/its_neighbor_index/ItsNeighborIndex.hpp b/sandboxes/its_neighbor_index/ItsNeighborIndex.hpp
new file mode 100644
index 0000000000..66f2b23485
--- /dev/null
+++ b/sandboxes/its_neighbor_index/ItsNeighborIndex.hpp
@@ -0,0 +1,19 @@
+#include <libslic3r/TriangleMesh.hpp>
+#include "libslic3r/MeshSplitImpl.hpp"
+
+namespace Slic3r {
+using FaceNeighborIndex = std::vector<std::array<size_t, 3>>;
+FaceNeighborIndex its_create_neighbors_index_1(const indexed_triangle_set &its);
+std::vector<Vec3i> its_create_neighbors_index_2(const indexed_triangle_set &its);
+std::vector<Vec3i> its_create_neighbors_index_3(const indexed_triangle_set &its);
+FaceNeighborIndex its_create_neighbors_index_4(const indexed_triangle_set &its);
+//FaceNeighborIndex its_create_neighbors_index_4(const indexed_triangle_set &its);
+std::vector<Vec3crd> its_create_neighbors_index_5(const indexed_triangle_set &its);
+std::vector<std::array<size_t, 3>> its_create_neighbors_index_6(const indexed_triangle_set &its);
+std::vector<std::array<size_t, 3>> its_create_neighbors_index_7(const indexed_triangle_set &its);
+FaceNeighborIndex its_create_neighbors_index_8(const indexed_triangle_set &its);
+std::vector<Vec3crd> its_create_neighbors_index_9(const indexed_triangle_set &its);
+std::vector<Vec3i> its_create_neighbors_index_10(const indexed_triangle_set &its);
+
+std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangle_set &its);
+}
diff --git a/sandboxes/its_neighbor_index/main.cpp b/sandboxes/its_neighbor_index/main.cpp
new file mode 100644
index 0000000000..2c89da3ac4
--- /dev/null
+++ b/sandboxes/its_neighbor_index/main.cpp
@@ -0,0 +1,273 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <tuple>
+#include <random>
+
+#include "ItsNeighborIndex.hpp"
+
+#include "libnest2d/tools/benchmark.h"
+#include "libnest2d/utils/metaloop.hpp"
+
+namespace Slic3r {
+
+enum { IndexCreation, Split };
+struct MeasureResult
+{
+    static constexpr const char * Names[] = {
+        "Index creation [s]",
+        "Split [s]"
+    };
+
+    double measurements[std::size(Names)] = {0.};
+};
+
+template<class IndexCreatorFn>
+static MeasureResult measure_index(const indexed_triangle_set &its, IndexCreatorFn fn)
+{
+    Benchmark b;
+
+    MeasureResult r;
+    for (int i = 0; i < 10; ++i) {
+        b.start();
+        ItsNeighborsWrapper itsn{its, fn(its)};
+        b.stop();
+
+        r.measurements[IndexCreation] += b.getElapsedSec();
+
+        b.start();
+        auto res = its_split(itsn);
+        b.stop();
+
+//        if (res.size() != 2 || res[0].indices.size() != res[1].indices.size() )
+//            std::cerr << "Something is wrong, split result invalid" << std::endl;
+
+        r.measurements[Split] += b.getElapsedSec();
+    }
+
+    r.measurements[IndexCreation] /= 10;
+    r.measurements[Split] /= 10;
+
+    return r;
+}
+
+const auto Seed = 0;// std::random_device{}();
+
+static indexed_triangle_set make_sphere_rnd(double radius, double detail)
+{
+    using namespace Slic3r;
+
+    auto sphere = its_make_sphere(radius, detail);
+
+    auto vfidx = create_vertex_faces_index(sphere);
+
+    const size_t vertexnum = sphere.vertices.size();
+    const size_t facenum   = sphere.indices.size();
+
+    std::mt19937 rng{Seed};
+    std::uniform_int_distribution<size_t> distv(sphere.vertices.size() / 2, sphere.vertices.size() - 1);
+    std::uniform_int_distribution<size_t> distf(sphere.indices.size() / 2, sphere.indices.size() - 1) ;
+
+    std::vector<bool> was(vertexnum / 2, false);
+
+    for (size_t i = 0; i < vertexnum / 2; ++i) {
+        size_t image = distv(rng);
+        if (was[image - vertexnum / 2]) continue;
+        was[image - vertexnum / 2] = true;
+
+        std::swap(sphere.vertices[i], sphere.vertices[image]);
+        for (size_t face_id : vfidx[i]) {
+            for (int &vi : sphere.indices[face_id])
+                if (vi == int(i)) vi = image;
+        }
+
+        for (size_t face_id : vfidx[image]) {
+            for (int &vi : sphere.indices[face_id])
+                if (vi == int(image)) vi = i;
+        }
+
+        std::swap(vfidx[i], vfidx[image]);
+    }
+
+    for (size_t i = 0; i < facenum / 2; ++i) {
+        size_t image = distf(rng);
+        std::swap(sphere.indices[i], sphere.indices[image]);
+    }
+
+    return sphere;
+}
+
+static indexed_triangle_set two_spheres(double detail)
+{
+    auto sphere1 = make_sphere_rnd(10., 2 * PI / detail), sphere2 = sphere1;
+
+    its_transform(sphere1, identity3f().translate(Vec3f{-5.f, 0.f, 0.f}));
+    its_transform(sphere2, identity3f().translate(Vec3f{5.f, 0.f, 0.f}));
+
+    its_merge(sphere1, sphere2);
+
+    return sphere1;
+}
+
+static indexed_triangle_set make_spheres(unsigned N, double detail)
+{
+    indexed_triangle_set ret, sphere = make_sphere_rnd(10., 2. * PI / detail);
+
+    for (unsigned i = 0u ; i < N; ++i)
+        its_merge(ret, sphere);
+
+    return ret;
+}
+
+constexpr double sq2 = std::sqrt(2.);
+
+static const std::pair<const std::string, indexed_triangle_set> ToMeasure[] = {
+//    {"two_spheres_1x", two_spheres(60.)},
+//    {"two_spheres_2x", two_spheres(120.)},
+//    {"two_spheres_4x", two_spheres(240.)},
+//    {"two_spheres_8x", two_spheres(480.)},
+//    {"two_spheres_16x", two_spheres(2 * 480.)},
+//    {"two_spheres_32x", two_spheres(2 * 2 * 480.)},
+
+    {"two_spheres_1x", two_spheres(60.)},
+    {"two_spheres_2x", two_spheres(sq2 * 60.)},
+    {"two_spheres_4x", two_spheres(2 * 60.)},
+    {"two_spheres_8x", two_spheres(sq2 * 2. * 60.)},
+    {"two_spheres_16x", two_spheres(4. * 60.)},
+    {"two_spheres_32x", two_spheres(sq2 * 4. * 60.)},
+    {"two_spheres_64x", two_spheres(8. * 60.)},
+    {"two_spheres_128x", two_spheres(sq2 * 8. * 60.)},
+    {"two_spheres_256x", two_spheres(16. * 60.)},
+    {"two_spheres_512x", two_spheres(sq2 * 16. * 60.)},
+
+    {"2_spheres", make_spheres(2, 60.)},
+    {"4_spheres", make_spheres(4, 60.)},
+    {"8_spheres", make_spheres(8, 60.)},
+    {"16_spheres", make_spheres(16,  60.)},
+    {"32_spheres", make_spheres(32, 60.)},
+    {"64_spheres", make_spheres(64, 60.)},
+    {"128_spheres", make_spheres(128, 60.)},
+    {"256_spheres", make_spheres(256, 60.)},
+    {"512_spheres", make_spheres(512, 60.)},
+    {"1024_spheres", make_spheres(1024, 60.)}
+};
+
+static const auto IndexFunctions = std::make_tuple(
+    std::make_pair("tamas's unordered_map based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_1); }),
+    std::make_pair("vojta std::sort based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_2); }),
+    std::make_pair("vojta tbb::parallel_sort based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_3); }),
+    std::make_pair("filip's vertex->face based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_5); }),
+    std::make_pair("vojta's vertex->face", [](const auto &its) { return measure_index(its, its_create_neighbors_index_9); }),
+    std::make_pair("vojta's vertex->face parallel", [](const auto &its) { return measure_index(its, its_create_neighbors_index_10); }),
+    std::make_pair("tamas's std::sort based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_6); }),
+    std::make_pair("tamas's tbb::parallel_sort based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_7); }),
+    std::make_pair("tamas's map based", [](const auto &its) { return measure_index(its, its_create_neighbors_index_8); }),
+    std::make_pair("TriangleMesh split", [](const auto &its) {
+
+        MeasureResult r;
+        for (int i = 0; i < 10; ++i) {
+            TriangleMesh m{its};
+            Benchmark b;
+
+            b.start();
+            m.repair(); // FIXME: this does more than just create neighborhood map
+            b.stop();
+            r.measurements[IndexCreation] += b.getElapsedSec();
+
+            b.start();
+            auto res = m.split();
+            b.stop();
+            r.measurements[Split] += b.getElapsedSec();
+
+//            if (res.size() != 2 || res[0]->size() != res[1]->size())
+//                std::cerr << "Something is wrong, split result invalid" << std::endl;
+        }
+        r.measurements[IndexCreation] /= 10;
+        r.measurements[Split] /= 10;
+
+        return r;
+    })
+
+//    std::make_pair("Vojta's vertex->face index", [](const auto &its){
+//        Benchmark b;
+//        b.start();
+//        auto index = create_vertex_faces_index(its);
+//        b.stop();
+
+//        if (index.size() != its.vertices.size())
+//            std::cerr << "Something went wrong!";
+
+//        return  MeasureResult{b.getElapsedSec(), 0., 0.};
+//    }),
+//    std::make_pair("Tamas's vertex->face index", [](const auto &its){
+//        Benchmark b;
+//        b.start();
+//        VertexFaceIndex index{its};
+//        b.stop();
+
+//        if (index.size() < its.vertices.size())
+//            std::cerr << "Something went wrong!";
+
+//        return  MeasureResult{b.getElapsedSec(), 0., 0.};
+//    })
+);
+
+static constexpr size_t IndexFuncNum = std::tuple_size_v<decltype (IndexFunctions)>;
+
+} // namespace Slic3r
+
+int main(const int argc, const char * argv[])
+{
+    using namespace Slic3r;
+
+    std::array<MeasureResult, IndexFuncNum> results[std::size(ToMeasure)];
+    std::array<std::string, IndexFuncNum> funcnames;
+
+    for (size_t i = 0; i < std::size(ToMeasure); ++i) {
+        auto &m = ToMeasure[i];
+        auto &name = m.first;
+        auto &mesh = m.second;
+
+//        its_write_obj(mesh, (std::string(name) + ".obj").c_str());
+
+        std::cout << "Mesh " << name << " has " << mesh.indices.size() << " faces and " << mesh.vertices.size() << " vertices." << std::endl;
+        libnest2d::opt::metaloop::apply([&mesh, i, &results, &funcnames](int N, auto &e) {
+            MeasureResult r = e.second(mesh);
+            funcnames[N] = e.first;
+            results[i][N] = r;
+        }, IndexFunctions);
+    }
+
+    std::string outfilename = "out.csv";
+    std::fstream outfile;
+    if (argc > 1) {
+        outfilename = argv[1];
+        outfile.open(outfilename, std::fstream::out);
+        std::cout << outfilename << " will be used" << std::endl;
+    }
+
+    std::ostream &out = outfile.is_open() ? outfile : std::cout;
+
+    for (size_t m = 0; m < std::size(MeasureResult::Names); ++m) {
+        out << MeasureResult::Names[m] << "\n";
+        out << std::endl;
+        out << "model;" ;
+        for (const std::string &funcname : funcnames) {
+            out << funcname << ";";
+        }
+
+        out << std::endl;
+
+        for (size_t i = 0; i < std::size(ToMeasure); ++i) {
+            const auto &result_row = results[i];
+            const std::string &name = ToMeasure[i].first;
+            out << name << ";";
+            for (auto &r : result_row)
+                out << r.measurements[m] << ";";
+
+            out << std::endl;
+        }
+    }
+
+    return 0;
+}
diff --git a/src/admesh/stl.h b/src/admesh/stl.h
index f7b8937ade..1ca75c95d4 100644
--- a/src/admesh/stl.h
+++ b/src/admesh/stl.h
@@ -149,9 +149,11 @@ 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
+
+    bool empty() const { return indices.empty() || vertices.empty(); }
 };
 
 extern bool stl_open(stl_file *stl, const char *file);
diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt
index 47e9ff5aed..2b5f3f9324 100644
--- a/src/libslic3r/CMakeLists.txt
+++ b/src/libslic3r/CMakeLists.txt
@@ -205,6 +205,7 @@ add_library(libslic3r STATIC
     TriangleMesh.hpp
     TriangleMeshSlicer.cpp
     TriangleMeshSlicer.hpp
+    MeshSplitImpl.hpp
     TriangulateWall.hpp
     TriangulateWall.cpp
     utils.cpp
@@ -266,8 +267,6 @@ add_library(libslic3r STATIC
     SLA/SupportPoint.hpp
     SLA/SupportPointGenerator.hpp
     SLA/SupportPointGenerator.cpp
-    SLA/Contour3D.hpp
-    SLA/Contour3D.cpp
     SLA/IndexedMesh.hpp
     SLA/IndexedMesh.cpp
     SLA/Clustering.hpp
diff --git a/src/libslic3r/Format/SL1.cpp b/src/libslic3r/Format/SL1.cpp
index 1809c3e4f3..f556b0ead7 100644
--- a/src/libslic3r/Format/SL1.cpp
+++ b/src/libslic3r/Format/SL1.cpp
@@ -296,7 +296,7 @@ void import_sla_archive(const std::string &zipfname, DynamicPrintConfig &out)
 void import_sla_archive(
     const std::string &      zipfname,
     Vec2i                    windowsize,
-    TriangleMesh &           out,
+    indexed_triangle_set &           out,
     DynamicPrintConfig &     profile,
     std::function<bool(int)> progr)
 {
@@ -316,7 +316,7 @@ void import_sla_archive(
         extract_slices_from_sla_archive(arch, rstp, progr);
 
     if (!slices.empty())
-        out = slices_to_triangle_mesh(slices, 0, slicp.layerh, slicp.initial_layerh);
+        out = slices_to_mesh(slices, 0, slicp.layerh, slicp.initial_layerh);
 }
 
 using ConfMap = std::map<std::string, std::string>;
diff --git a/src/libslic3r/Format/SL1.hpp b/src/libslic3r/Format/SL1.hpp
index ab731ff841..a3e6ce26cb 100644
--- a/src/libslic3r/Format/SL1.hpp
+++ b/src/libslic3r/Format/SL1.hpp
@@ -14,7 +14,7 @@ class SL1Archive: public SLAPrinter {
 protected:
     uqptr<sla::RasterBase> create_raster() const override;
     sla::RasterEncoder get_encoder() const override;
-    
+
 public:
     
     SL1Archive() = default;
@@ -43,14 +43,14 @@ void import_sla_archive(const std::string &zipfname, DynamicPrintConfig &out);
 void import_sla_archive(
     const std::string &      zipfname,
     Vec2i                    windowsize,
-    TriangleMesh &           out,
+    indexed_triangle_set &   out,
     DynamicPrintConfig &     profile,
     std::function<bool(int)> progr = [](int) { return true; });
 
 inline void import_sla_archive(
     const std::string &      zipfname,
     Vec2i                    windowsize,
-    TriangleMesh &           out,
+    indexed_triangle_set &   out,
     std::function<bool(int)> progr = [](int) { return true; })
 {
     DynamicPrintConfig profile;
diff --git a/src/libslic3r/MTUtils.hpp b/src/libslic3r/MTUtils.hpp
index 9e77aa90af..e60918fab6 100644
--- a/src/libslic3r/MTUtils.hpp
+++ b/src/libslic3r/MTUtils.hpp
@@ -74,29 +74,6 @@ public:
     }
 };
 
-/// A very simple range concept implementation with iterator-like objects.
-template<class It> class Range
-{
-    It from, to;
-
-public:
-    // The class is ready for range based for loops.
-    It begin() const { return from; }
-    It end() const { return to; }
-
-    // The iterator type can be obtained this way.
-    using Type = It;
-
-    Range() = default;
-    Range(It &&b, It &&e)
-        : from(std::forward<It>(b)), to(std::forward<It>(e))
-    {}
-
-    // Some useful container-like methods...
-    inline size_t size() const { return end() - begin(); }
-    inline bool   empty() const { return size() == 0; }
-};
-
 template<class C> bool all_of(const C &container)
 {
     return std::all_of(container.begin(),
diff --git a/src/libslic3r/MeshSplitImpl.hpp b/src/libslic3r/MeshSplitImpl.hpp
new file mode 100644
index 0000000000..322fc82f88
--- /dev/null
+++ b/src/libslic3r/MeshSplitImpl.hpp
@@ -0,0 +1,226 @@
+#ifndef MESHSPLITIMPL_HPP
+#define MESHSPLITIMPL_HPP
+
+#include "TriangleMesh.hpp"
+#include "libnest2d/tools/benchmark.h"
+#include "Execution/ExecutionTBB.hpp"
+
+namespace Slic3r {
+
+template<class ExPolicy>
+std::vector<Vec3i> create_neighbors_index(ExPolicy &&ex, const indexed_triangle_set &its);
+
+namespace meshsplit_detail {
+
+template<class Its, class Enable = void> struct ItsWithNeighborsIndex_ {
+    using Index = typename Its::Index;
+    static const indexed_triangle_set &get_its(const Its &m) { return m.get_its();}
+    static const Index &get_index(const Its &m) { return m.get_index(); }
+};
+
+// Define a default neighbors index for indexed_triangle_set
+template<> struct ItsWithNeighborsIndex_<indexed_triangle_set> {
+    using Index = std::vector<Vec3i>;
+    static const indexed_triangle_set &get_its(const indexed_triangle_set &its) noexcept { return its; }
+    static Index get_index(const indexed_triangle_set &its) noexcept
+    {
+        return create_neighbors_index(ex_tbb, its);
+    }
+};
+
+// Visit all unvisited neighboring facets that are reachable from the first unvisited facet,
+// and return them.
+template<class NeighborIndex>
+std::vector<size_t> its_find_unvisited_neighbors(
+    const indexed_triangle_set &its,
+    const NeighborIndex &       neighbor_index,
+    std::vector<char> &         visited)
+{
+    using stack_el = size_t;
+
+    auto facestack = reserve_vector<stack_el>(its.indices.size());
+    auto push = [&facestack] (const stack_el &s) { facestack.emplace_back(s); };
+    auto pop  = [&facestack] () -> stack_el {
+        stack_el ret = facestack.back();
+        facestack.pop_back();
+        return ret;
+    };
+
+    // find the next unvisited facet and push the index
+    auto facet = std::find(visited.begin(), visited.end(), false);
+    std::vector<size_t> ret;
+
+    if (facet != visited.end()) {
+        ret.reserve(its.indices.size());
+        auto idx = size_t(facet - visited.begin());
+        push(idx);
+        ret.emplace_back(idx);
+        visited[idx] = true;
+    }
+
+    while (!facestack.empty()) {
+        size_t facet_idx = pop();
+        const auto &neighbors = neighbor_index[facet_idx];
+        for (auto neighbor_idx : neighbors) {
+            if (size_t(neighbor_idx) < visited.size() && !visited[size_t(neighbor_idx)]) {
+                visited[size_t(neighbor_idx)] = true;
+                push(stack_el(neighbor_idx));
+                ret.emplace_back(size_t(neighbor_idx));
+            }
+        }
+    }
+
+    return ret;
+}
+
+} // namespace meshsplit_detail
+
+template<class IndexT> struct ItsNeighborsWrapper
+{
+    using Index = IndexT;
+    const indexed_triangle_set *its;
+    IndexT index;
+
+    ItsNeighborsWrapper(const indexed_triangle_set &m, IndexT &&idx)
+        : its{&m}, index{std::move(idx)}
+    {}
+
+    const auto& get_its() const noexcept { return *its; }
+    const auto& get_index() const noexcept { return index; }
+};
+
+// Splits a mesh into multiple meshes when possible.
+template<class Its, class OutputIt>
+void its_split(const Its &m, OutputIt out_it)
+{
+    using namespace meshsplit_detail;
+
+    const indexed_triangle_set &its = ItsWithNeighborsIndex_<Its>::get_its(m);
+
+    std::vector<char> visited(its.indices.size(), false);
+
+    struct VertexConv {
+        size_t part_id      = std::numeric_limits<size_t>::max();
+        size_t vertex_image;
+    };
+    std::vector<VertexConv> vidx_conv(its.vertices.size());
+
+    const auto& neighbor_index = ItsWithNeighborsIndex_<Its>::get_index(m);
+
+    for (size_t part_id = 0;; ++part_id) {
+        std::vector<size_t> facets =
+            its_find_unvisited_neighbors(its, neighbor_index, visited);
+
+        if (facets.empty())
+            break;
+
+        // Create a new mesh for the part that was just split off.
+        indexed_triangle_set mesh;
+        mesh.indices.reserve(facets.size());
+        mesh.vertices.reserve(std::min(facets.size() * 3, its.vertices.size()));
+
+        // Assign the facets to the new mesh.
+        for (size_t face_id : facets) {
+            const auto &face = its.indices[face_id];
+            Vec3i       new_face;
+            for (size_t v = 0; v < 3; ++v) {
+                auto vi = face(v);
+
+                if (vidx_conv[vi].part_id != part_id) {
+                    vidx_conv[vi] = {part_id, mesh.vertices.size()};
+                    mesh.vertices.emplace_back(its.vertices[size_t(vi)]);
+                }
+
+                new_face(v) = vidx_conv[vi].vertex_image;
+            }
+
+            mesh.indices.emplace_back(new_face);
+        }
+
+        out_it = std::move(mesh);
+    }
+}
+
+template<class Its>
+std::vector<indexed_triangle_set> its_split(const Its &its)
+{
+    auto ret = reserve_vector<indexed_triangle_set>(3);
+    its_split(its, std::back_inserter(ret));
+
+    return ret;
+}
+
+template<class Its> bool its_is_splittable(const Its &m)
+{
+    using namespace meshsplit_detail;
+    const indexed_triangle_set &its = ItsWithNeighborsIndex_<Its>::get_its(m);
+    const auto& neighbor_index = ItsWithNeighborsIndex_<Its>::get_index(m);
+
+    std::vector<char> visited(its.indices.size(), false);
+    auto faces = its_find_unvisited_neighbors(its, neighbor_index, visited);
+
+    return !faces.empty();
+}
+
+inline int get_vertex_index(size_t vertex_index, const stl_triangle_vertex_indices &triangle_indices) {
+    if (int(vertex_index) == triangle_indices[0]) return 0;
+    if (int(vertex_index) == triangle_indices[1]) return 1;
+    if (int(vertex_index) == triangle_indices[2]) return 2;
+    return -1;
+}
+
+inline Vec2crd get_edge_indices(int edge_index, const stl_triangle_vertex_indices &triangle_indices)
+{
+    int next_edge_index = (edge_index == 2) ? 0 : edge_index + 1;
+    int vi0             = triangle_indices[edge_index];
+    int vi1             = triangle_indices[next_edge_index];
+    return Vec2crd(vi0, vi1);
+}
+
+template<class ExPolicy>
+std::vector<Vec3i> create_neighbors_index(ExPolicy &&ex, const indexed_triangle_set &its)
+{
+    const std::vector<stl_triangle_vertex_indices> &indices = its.indices;
+    size_t vertices_size = its.vertices.size();
+
+    if (indices.empty() || vertices_size == 0) return {};
+
+    auto               vertex_triangles = VertexFaceIndex{its};
+    static constexpr int no_value         = -1;
+    std::vector<Vec3i> neighbors(indices.size(),
+                                 Vec3i(no_value, no_value, no_value));
+
+    //for (const stl_triangle_vertex_indices& triangle_indices : indices) {
+    execution::for_each(ex, size_t(0), indices.size(),
+        [&neighbors, &indices, &vertex_triangles] (size_t index)
+        {
+            Vec3i& neighbor = neighbors[index];
+            const stl_triangle_vertex_indices & triangle_indices = indices[index];
+            for (int edge_index = 0; edge_index < 3; ++edge_index) {
+                // check if done
+                int& neighbor_edge = neighbor[edge_index];
+                if (neighbor_edge != no_value) continue;
+                Vec2crd edge_indices = get_edge_indices(edge_index, triangle_indices);
+                // IMPROVE: use same vector for 2 sides of triangle
+                const auto &faces_range = vertex_triangles[edge_indices[0]];
+                for (const size_t &face : faces_range) {
+                    if (face <= index) continue;
+                    const stl_triangle_vertex_indices &face_indices = indices[face];
+                    int vertex_index = get_vertex_index(edge_indices[1], face_indices);
+                    // NOT Contain second vertex?
+                    if (vertex_index < 0) continue;
+                    // Has NOT oposit direction?
+                    if (edge_indices[0] != face_indices[(vertex_index + 1) % 3]) continue;
+                    neighbor_edge = face;
+                    neighbors[face][vertex_index] = index;
+                    break;
+                }
+            }
+        }, execution::max_concurrency(ex));
+
+    return neighbors;
+}
+
+} // namespace Slic3r
+
+#endif // MESHSPLITIMPL_HPP
diff --git a/src/libslic3r/OpenVDBUtils.cpp b/src/libslic3r/OpenVDBUtils.cpp
index e09fceed7a..72c7668a45 100644
--- a/src/libslic3r/OpenVDBUtils.cpp
+++ b/src/libslic3r/OpenVDBUtils.cpp
@@ -21,11 +21,11 @@ namespace Slic3r {
 
 class TriangleMeshDataAdapter {
 public:
-    const TriangleMesh &mesh;
+    const indexed_triangle_set &its;
     float voxel_scale;
 
-    size_t polygonCount() const { return mesh.its.indices.size(); }
-    size_t pointCount() const   { return mesh.its.vertices.size(); }
+    size_t polygonCount() const { return its.indices.size(); }
+    size_t pointCount() const   { return its.vertices.size(); }
     size_t vertexCount(size_t) const { return 3; }
 
     // Return position pos in local grid index space for polygon n and vertex v
@@ -33,19 +33,19 @@ public:
     // And the voxel count per unit volume can be affected this way.
     void getIndexSpacePoint(size_t n, size_t v, openvdb::Vec3d& pos) const
     {
-        auto vidx = size_t(mesh.its.indices[n](Eigen::Index(v)));
-        Slic3r::Vec3d p = mesh.its.vertices[vidx].cast<double>() * voxel_scale;
+        auto vidx = size_t(its.indices[n](Eigen::Index(v)));
+        Slic3r::Vec3d p = its.vertices[vidx].cast<double>() * voxel_scale;
         pos = {p.x(), p.y(), p.z()};
     }
 
-    TriangleMeshDataAdapter(const TriangleMesh &m, float voxel_sc = 1.f)
-        : mesh{m}, voxel_scale{voxel_sc} {};
+    TriangleMeshDataAdapter(const indexed_triangle_set &m, float voxel_sc = 1.f)
+        : its{m}, voxel_scale{voxel_sc} {};
 };
 
 // TODO: Do I need to call initialize? Seems to work without it as well but the
 // docs say it should be called ones. It does a mutex lock-unlock sequence all
 // even if was called previously.
-openvdb::FloatGrid::Ptr mesh_to_grid(const TriangleMesh &            mesh,
+openvdb::FloatGrid::Ptr mesh_to_grid(const indexed_triangle_set &    mesh,
                                      const openvdb::math::Transform &tr,
                                      float voxel_scale,
                                      float exteriorBandWidth,
@@ -54,22 +54,17 @@ openvdb::FloatGrid::Ptr mesh_to_grid(const TriangleMesh &            mesh,
 {
     openvdb::initialize();
 
-    TriangleMeshPtrs meshparts_raw = mesh.split();
-    auto meshparts = reserve_vector<std::unique_ptr<TriangleMesh>>(meshparts_raw.size());
-    for (auto *p : meshparts_raw)
-        meshparts.emplace_back(p);
+    std::vector<indexed_triangle_set> meshparts = its_split(mesh);
 
-    auto it = std::remove_if(meshparts.begin(), meshparts.end(), [](auto &m) {
-         m->require_shared_vertices();
-         return m->volume() < EPSILON;
-     });
+    auto it = std::remove_if(meshparts.begin(), meshparts.end(),
+                             [](auto &m) { return its_volume(m) < EPSILON; });
 
     meshparts.erase(it, meshparts.end());
 
     openvdb::FloatGrid::Ptr grid;
     for (auto &m : meshparts) {
         auto subgrid = openvdb::tools::meshToVolume<openvdb::FloatGrid>(
-            TriangleMeshDataAdapter{*m, voxel_scale}, tr, exteriorBandWidth,
+            TriangleMeshDataAdapter{m, voxel_scale}, tr, exteriorBandWidth,
             interiorBandWidth, flags);
 
         if (grid && subgrid) openvdb::tools::csgUnion(*grid, *subgrid);
@@ -91,11 +86,10 @@ openvdb::FloatGrid::Ptr mesh_to_grid(const TriangleMesh &            mesh,
     return grid;
 }
 
-template<class Grid>
-sla::Contour3D _volumeToMesh(const Grid &grid,
-                             double      isovalue,
-                             double      adaptivity,
-                             bool        relaxDisorientedTriangles)
+indexed_triangle_set grid_to_mesh(const openvdb::FloatGrid &grid,
+                          double                    isovalue,
+                          double                    adaptivity,
+                          bool                      relaxDisorientedTriangles)
 {
     openvdb::initialize();
 
@@ -111,36 +105,20 @@ sla::Contour3D _volumeToMesh(const Grid &grid,
         scale = grid.template metaValue<float>("voxel_scale");
     }  catch (...) { }
 
-    sla::Contour3D ret;
-    ret.points.reserve(points.size());
-    ret.faces3.reserve(triangles.size());
-    ret.faces4.reserve(quads.size());
+    indexed_triangle_set ret;
+    ret.vertices.reserve(points.size());
+    ret.indices.reserve(triangles.size() + quads.size() * 2);
 
-    for (auto &v : points) ret.points.emplace_back(to_vec3d(v) / scale);
-    for (auto &v : triangles) ret.faces3.emplace_back(to_vec3i(v));
-    for (auto &v : quads) ret.faces4.emplace_back(to_vec4i(v));
+    for (auto &v : points) ret.vertices.emplace_back(to_vec3f(v) / scale);
+    for (auto &v : triangles) ret.indices.emplace_back(to_vec3i(v));
+    for (auto &quad : quads) {
+        ret.indices.emplace_back(quad(0), quad(1), quad(2));
+        ret.indices.emplace_back(quad(2), quad(3), quad(0));
+    }
 
     return ret;
 }
 
-TriangleMesh grid_to_mesh(const openvdb::FloatGrid &grid,
-                          double                    isovalue,
-                          double                    adaptivity,
-                          bool                      relaxDisorientedTriangles)
-{
-    return to_triangle_mesh(
-        _volumeToMesh(grid, isovalue, adaptivity, relaxDisorientedTriangles));
-}
-
-sla::Contour3D grid_to_contour3d(const openvdb::FloatGrid &grid,
-                                 double                    isovalue,
-                                 double                    adaptivity,
-                                 bool relaxDisorientedTriangles)
-{
-    return _volumeToMesh(grid, isovalue, adaptivity,
-                         relaxDisorientedTriangles);
-}
-
 openvdb::FloatGrid::Ptr redistance_grid(const openvdb::FloatGrid &grid,
                                         double                    iso,
                                         double                    er,
diff --git a/src/libslic3r/OpenVDBUtils.hpp b/src/libslic3r/OpenVDBUtils.hpp
index 5df816abb1..490b3e2228 100644
--- a/src/libslic3r/OpenVDBUtils.hpp
+++ b/src/libslic3r/OpenVDBUtils.hpp
@@ -2,7 +2,6 @@
 #define OPENVDBUTILS_HPP
 
 #include <libslic3r/TriangleMesh.hpp>
-#include <libslic3r/SLA/Contour3D.hpp>
 
 #ifdef _MSC_VER
 // Suppress warning C4146 in include/gmp.h(2177,31): unary minus operator applied to unsigned type, result still unsigned 
@@ -19,7 +18,6 @@ namespace Slic3r {
 inline Vec3f to_vec3f(const openvdb::Vec3s &v) { return Vec3f{v.x(), v.y(), v.z()}; }
 inline Vec3d to_vec3d(const openvdb::Vec3s &v) { return to_vec3f(v).cast<double>(); }
 inline Vec3i to_vec3i(const openvdb::Vec3I &v) { return Vec3i{int(v[0]), int(v[1]), int(v[2])}; }
-inline Vec4i to_vec4i(const openvdb::Vec4I &v) { return Vec4i{int(v[0]), int(v[1]), int(v[2]), int(v[3])}; }
 
 // Here voxel_scale defines the scaling of voxels which affects the voxel count.
 // 1.0 value means a voxel for every unit cube. 2 means the model is scaled to
@@ -28,22 +26,17 @@ inline Vec4i to_vec4i(const openvdb::Vec4I &v) { return Vec4i{int(v[0]), int(v[1
 // achievable through the Transform parameter. (TODO: or is it?)
 // The resulting grid will contain the voxel_scale in its metadata under the
 // "voxel_scale" key to be used in grid_to_mesh function.
-openvdb::FloatGrid::Ptr mesh_to_grid(const TriangleMesh &            mesh,
+openvdb::FloatGrid::Ptr mesh_to_grid(const indexed_triangle_set &    mesh,
                                      const openvdb::math::Transform &tr = {},
-                                     float voxel_scale = 1.f,
+                                     float voxel_scale                  = 1.f,
                                      float exteriorBandWidth = 3.0f,
                                      float interiorBandWidth = 3.0f,
                                      int   flags             = 0);
 
-sla::Contour3D grid_to_contour3d(const openvdb::FloatGrid &grid,
-                                 double                    isovalue,
-                                 double                    adaptivity,
-                                 bool relaxDisorientedTriangles = true);
-
-TriangleMesh grid_to_mesh(const openvdb::FloatGrid &grid,
-                          double                    isovalue   = 0.0,
-                          double                    adaptivity = 0.0,
-                          bool relaxDisorientedTriangles       = true);
+indexed_triangle_set grid_to_mesh(const openvdb::FloatGrid &grid,
+                                  double                    isovalue   = 0.0,
+                                  double                    adaptivity = 0.0,
+                                  bool relaxDisorientedTriangles = true);
 
 openvdb::FloatGrid::Ptr redistance_grid(const openvdb::FloatGrid &grid,
                                         double                    iso,
diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp
index 248d4129f4..20726270d8 100644
--- a/src/libslic3r/Point.hpp
+++ b/src/libslic3r/Point.hpp
@@ -21,6 +21,12 @@ class MultiPoint;
 class Point;
 using Vector = Point;
 
+// Base template for eigen derived vectors
+template<int N, int M, class T>
+using Mat = Eigen::Matrix<T, N, M, Eigen::DontAlign, N, M>;
+
+template<int N, class T> using Vec = Mat<N, 1, T>;
+
 // Eigen types, to replace the Slic3r's own types in the future.
 // Vector types with a fixed point coordinate base type.
 using Vec2crd = Eigen::Matrix<coord_t,  2, 1, Eigen::DontAlign>;
@@ -54,11 +60,19 @@ using Matrix3d       = Eigen::Matrix<double, 3, 3, Eigen::DontAlign>;
 using Matrix4f       = Eigen::Matrix<float,  4, 4, Eigen::DontAlign>;
 using Matrix4d       = Eigen::Matrix<double, 4, 4, Eigen::DontAlign>;
 
+template<int N, class T>
+using Transform = Eigen::Transform<float, N, Eigen::Affine, Eigen::DontAlign>;
+
 using Transform2f    = Eigen::Transform<float,  2, Eigen::Affine, Eigen::DontAlign>;
 using Transform2d    = Eigen::Transform<double, 2, Eigen::Affine, Eigen::DontAlign>;
 using Transform3f    = Eigen::Transform<float,  3, Eigen::Affine, Eigen::DontAlign>;
 using Transform3d    = Eigen::Transform<double, 3, Eigen::Affine, Eigen::DontAlign>;
 
+// I don't know why Eigen::Transform::Identity() return a const object...
+template<int N, class T> Transform<N, T> identity() { return Transform<N, T>::Identity(); }
+inline const auto &identity3f = identity<3, float>;
+inline const auto &identity3d = identity<3, double>;
+
 inline bool operator<(const Vec2d &lhs, const Vec2d &rhs) { return lhs.x() < rhs.x() || (lhs.x() == rhs.x() && lhs.y() < rhs.y()); }
 
 template<int Options>
@@ -488,4 +502,19 @@ namespace cereal {
 	template<class Archive> void save(Archive& archive, Slic3r::Matrix2f &m) { archive.saveBinary((char*)m.data(), sizeof(float) * 4); }
 }
 
+// To be able to use Vec<> and Mat<> in range based for loops:
+namespace Eigen {
+template<class T, int N, int M>
+T* begin(Slic3r::Mat<N, M, T> &mat) { return mat.data(); }
+
+template<class T, int N, int M>
+T* end(Slic3r::Mat<N, M, T> &mat) { return mat.data() + N * M; }
+
+template<class T, int N, int M>
+const T* begin(const Slic3r::Mat<N, M, T> &mat) { return mat.data(); }
+
+template<class T, int N, int M>
+const T* end(const Slic3r::Mat<N, M, T> &mat) { return mat.data() + N * M; }
+} // namespace Eigen
+
 #endif
diff --git a/src/libslic3r/SLA/Contour3D.cpp b/src/libslic3r/SLA/Contour3D.cpp
deleted file mode 100644
index 96d10af208..0000000000
--- a/src/libslic3r/SLA/Contour3D.cpp
+++ /dev/null
@@ -1,149 +0,0 @@
-#include <libslic3r/SLA/Contour3D.hpp>
-#include <libslic3r/SLA/IndexedMesh.hpp>
-
-#include <libslic3r/Format/objparser.hpp>
-
-namespace Slic3r { namespace sla {
-
-Contour3D::Contour3D(const TriangleMesh &trmesh)
-{
-    points.reserve(trmesh.its.vertices.size());
-    faces3.reserve(trmesh.its.indices.size());
-    
-    for (auto &v : trmesh.its.vertices)
-        points.emplace_back(v.cast<double>());
-    
-    std::copy(trmesh.its.indices.begin(), trmesh.its.indices.end(),
-              std::back_inserter(faces3));
-}
-
-Contour3D::Contour3D(TriangleMesh &&trmesh)
-{
-    points.reserve(trmesh.its.vertices.size());
-    
-    for (auto &v : trmesh.its.vertices)
-        points.emplace_back(v.cast<double>());
-    
-    faces3.swap(trmesh.its.indices);
-}
-
-Contour3D::Contour3D(const IndexedMesh &emesh) {
-    points.reserve(emesh.vertices().size());
-    faces3.reserve(emesh.indices().size());
-    
-    for (const Vec3f& vert : emesh.vertices())
-        points.emplace_back(vert.cast<double>());
-    
-    for (const auto& ind : emesh.indices())
-        faces3.emplace_back(ind);
-}
-
-Contour3D &Contour3D::merge(const Contour3D &ctr)
-{
-    auto N = coord_t(points.size());
-    auto N_f3 = faces3.size();
-    auto N_f4 = faces4.size();
-    
-    points.insert(points.end(), ctr.points.begin(), ctr.points.end());
-    faces3.insert(faces3.end(), ctr.faces3.begin(), ctr.faces3.end());
-    faces4.insert(faces4.end(), ctr.faces4.begin(), ctr.faces4.end());
-    
-    for(size_t n = N_f3; n < faces3.size(); n++) {
-        auto& idx = faces3[n]; idx.x() += N; idx.y() += N; idx.z() += N;
-    }
-    
-    for(size_t n = N_f4; n < faces4.size(); n++) {
-        auto& idx = faces4[n]; for (int k = 0; k < 4; k++) idx(k) += N;
-    }        
-    
-    return *this;
-}
-
-Contour3D &Contour3D::merge(const Pointf3s &triangles)
-{
-    const size_t offs = points.size();
-    points.insert(points.end(), triangles.begin(), triangles.end());
-    faces3.reserve(faces3.size() + points.size() / 3);
-    
-    for(int i = int(offs); i < int(points.size()); i += 3)
-        faces3.emplace_back(i, i + 1, i + 2);
-    
-    return *this;
-}
-
-void Contour3D::to_obj(std::ostream &stream)
-{
-    for(auto& p : points)
-        stream << "v " << p.transpose() << "\n";
-    
-    for(auto& f : faces3) 
-        stream << "f " << (f + Vec3i(1, 1, 1)).transpose() << "\n";
-    
-    for(auto& f : faces4)
-        stream << "f " << (f + Vec4i(1, 1, 1, 1)).transpose() << "\n";
-}
-
-void Contour3D::from_obj(std::istream &stream)
-{
-    ObjParser::ObjData data;
-    ObjParser::objparse(stream, data);
-    
-    points.reserve(data.coordinates.size() / 4 + 1);
-    auto &coords = data.coordinates;
-    for (size_t i = 0; i < coords.size(); i += 4)
-        points.emplace_back(coords[i], coords[i + 1], coords[i + 2]);
-    
-    Vec3i triangle;
-    Vec4i quad;
-    size_t v = 0;
-    while(v < data.vertices.size()) {
-        size_t N = 0;
-        size_t i = v;
-        while (data.vertices[v++].coordIdx != -1) ++N;
-        
-        std::function<void(int, int)> setfn;
-        if (N < 3 || N > 4) continue;
-        else if (N == 3) setfn = [&triangle](int k, int f) { triangle(k) = f; };
-        else setfn = [&quad](int k, int f) { quad(k) = f; };
-        
-        for (size_t j = 0; j < N; ++j)
-            setfn(int(j), data.vertices[i + j].coordIdx);
-    }
-}
-
-TriangleMesh to_triangle_mesh(const Contour3D &ctour) {
-    if (ctour.faces4.empty()) return {ctour.points, ctour.faces3};
-    
-    std::vector<Vec3i> triangles;
-    
-    triangles.reserve(ctour.faces3.size() + 2 * ctour.faces4.size());
-    std::copy(ctour.faces3.begin(), ctour.faces3.end(),
-              std::back_inserter(triangles));
-    
-    for (auto &quad : ctour.faces4) {
-        triangles.emplace_back(quad(0), quad(1), quad(2));
-        triangles.emplace_back(quad(2), quad(3), quad(0));
-    }
-    
-    return {ctour.points, std::move(triangles)};
-}
-
-TriangleMesh to_triangle_mesh(Contour3D &&ctour) {
-    if (ctour.faces4.empty())
-        return {std::move(ctour.points), std::move(ctour.faces3)};
-    
-    std::vector<Vec3i> triangles;
-    
-    triangles.reserve(ctour.faces3.size() + 2 * ctour.faces4.size());
-    std::copy(ctour.faces3.begin(), ctour.faces3.end(),
-              std::back_inserter(triangles));
-    
-    for (auto &quad : ctour.faces4) {
-        triangles.emplace_back(quad(0), quad(1), quad(2));
-        triangles.emplace_back(quad(2), quad(3), quad(0));
-    }
-    
-    return {std::move(ctour.points), std::move(triangles)};
-}
-
-}} // namespace Slic3r::sla
diff --git a/src/libslic3r/SLA/Contour3D.hpp b/src/libslic3r/SLA/Contour3D.hpp
deleted file mode 100644
index 3380cd6ab0..0000000000
--- a/src/libslic3r/SLA/Contour3D.hpp
+++ /dev/null
@@ -1,48 +0,0 @@
-#ifndef SLA_CONTOUR3D_HPP
-#define SLA_CONTOUR3D_HPP
-
-#include <libslic3r/TriangleMesh.hpp>
-
-namespace Slic3r {
-
-// Used for quads (TODO: remove this, and convert quads to triangles in OpenVDBUtils)
-using Vec4i = Eigen::Matrix<int, 4, 1, Eigen::DontAlign>;
-
-namespace sla {
-
-class IndexedMesh;
-
-/// Dumb vertex mesh consisting of triangles (or) quads. Capable of merging with
-/// other meshes of this type and converting to and from other mesh formats.
-struct Contour3D {
-    std::vector<Vec3d> points;
-    std::vector<Vec3i> faces3;
-    std::vector<Vec4i> faces4;
-    
-    Contour3D() = default;
-    Contour3D(const TriangleMesh &trmesh);
-    Contour3D(TriangleMesh &&trmesh);
-    Contour3D(const IndexedMesh  &emesh);
-    
-    Contour3D& merge(const Contour3D& ctr);
-    Contour3D& merge(const Pointf3s& triangles);
-    
-    // Write the index triangle structure to OBJ file for debugging purposes.
-    void to_obj(std::ostream& stream);
-    void from_obj(std::istream &stream);
-
-    inline bool empty() const
-    {
-        return points.empty() || (faces4.empty() && faces3.empty());
-    }
-};
-
-/// Mesh from an existing contour.
-TriangleMesh to_triangle_mesh(const Contour3D& ctour);
-
-/// Mesh from an evaporating 3D contour
-TriangleMesh to_triangle_mesh(Contour3D&& ctour);
-
-}} // namespace Slic3r::sla
-
-#endif // CONTOUR3D_HPP
diff --git a/src/libslic3r/SLA/Hollowing.cpp b/src/libslic3r/SLA/Hollowing.cpp
index 32a8810375..770e52a359 100644
--- a/src/libslic3r/SLA/Hollowing.cpp
+++ b/src/libslic3r/SLA/Hollowing.cpp
@@ -22,14 +22,8 @@
 namespace Slic3r {
 namespace sla {
 
-template<class S, class = FloatingOnly<S>>
-inline void _scale(S s, TriangleMesh &m) { m.scale(float(s)); }
-
-template<class S, class = FloatingOnly<S>>
-inline void _scale(S s, Contour3D &m) { for (auto &p : m.points) p *= s; }
-
 struct Interior {
-    TriangleMesh mesh;
+    indexed_triangle_set mesh;
     openvdb::FloatGrid::Ptr gridptr;
     mutable std::optional<openvdb::FloatGrid::ConstAccessor> accessor;
 
@@ -53,12 +47,12 @@ void InteriorDeleter::operator()(Interior *p)
     delete p;
 }
 
-TriangleMesh &get_mesh(Interior &interior)
+indexed_triangle_set &get_mesh(Interior &interior)
 {
     return interior.mesh;
 }
 
-const TriangleMesh &get_mesh(const Interior &interior)
+const indexed_triangle_set &get_mesh(const Interior &interior)
 {
     return interior.mesh;
 }
@@ -77,7 +71,7 @@ static InteriorPtr generate_interior_verbose(const TriangleMesh & mesh,
     if (ctl.stopcondition()) return {};
     else ctl.statuscb(0, L("Hollowing"));
 
-    auto gridptr = mesh_to_grid(mesh, {}, voxel_scale, out_range, in_range);
+    auto gridptr = mesh_to_grid(mesh.its, {}, voxel_scale, out_range, in_range);
 
     assert(gridptr);
 
@@ -136,32 +130,28 @@ InteriorPtr generate_interior(const TriangleMesh &   mesh,
 
     if (interior && !interior->mesh.empty()) {
 
-        // This flips the normals to be outward facing...
-        interior->mesh.require_shared_vertices();
-        indexed_triangle_set its = std::move(interior->mesh.its);
+        // flip normals back...
+        swap_normals(interior->mesh);
+        Slic3r::simplify_mesh(interior->mesh);
 
-        Slic3r::simplify_mesh(its);
+        its_compactify_vertices(interior->mesh);
+        its_merge_vertices(interior->mesh);
 
         // flip normals back...
-        for (stl_triangle_vertex_indices &ind : its.indices)
-            std::swap(ind(0), ind(2));
-
-        interior->mesh = Slic3r::TriangleMesh{its};
-        interior->mesh.repaired = true;
-        interior->mesh.require_shared_vertices();
+        swap_normals(interior->mesh);
     }
 
     return interior;
 }
 
-Contour3D DrainHole::to_mesh() const
+indexed_triangle_set DrainHole::to_mesh() const
 {
     auto r = double(radius);
     auto h = double(height);
-    sla::Contour3D hole = sla::cylinder(r, h, steps);
-    Eigen::Quaterniond q;
-    q.setFromTwoVectors(Vec3d{0., 0., 1.}, normal.cast<double>());
-    for(auto& p : hole.points) p = q * p + pos.cast<double>();
+    indexed_triangle_set hole = sla::cylinder(r, h, steps);
+    Eigen::Quaternionf q;
+    q.setFromTwoVectors(Vec3f{0.f, 0.f, 1.f}, normal);
+    for(auto& p : hole.vertices) p = q * p + pos;
     
     return hole;
 }
@@ -292,7 +282,7 @@ void cut_drainholes(std::vector<ExPolygons> & obj_slices,
 {
     TriangleMesh mesh;
     for (const sla::DrainHole &holept : holes)
-        mesh.merge(sla::to_triangle_mesh(holept.to_mesh()));
+        mesh.merge(TriangleMesh{holept.to_mesh()});
     
     if (mesh.empty()) return;
     
@@ -325,7 +315,7 @@ void hollow_mesh(TriangleMesh &mesh, const Interior &interior, int flags)
     if (flags & hfRemoveInsideTriangles && interior.gridptr)
         remove_inside_triangles(mesh, interior);
 
-    mesh.merge(interior.mesh);
+    mesh.merge(TriangleMesh{interior.mesh});
     mesh.require_shared_vertices();
 }
 
diff --git a/src/libslic3r/SLA/Hollowing.hpp b/src/libslic3r/SLA/Hollowing.hpp
index d9a77cd35e..b57513fe72 100644
--- a/src/libslic3r/SLA/Hollowing.hpp
+++ b/src/libslic3r/SLA/Hollowing.hpp
@@ -2,13 +2,11 @@
 #define SLA_HOLLOWING_HPP
 
 #include <memory>
-#include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/TriangleMesh.hpp>
 #include <libslic3r/SLA/JobController.hpp>
 
 namespace Slic3r {
 
-class TriangleMesh;
-
 namespace sla {
 
 struct HollowingConfig
@@ -27,8 +25,8 @@ struct Interior;
 struct InteriorDeleter { void operator()(Interior *p); };
 using  InteriorPtr = std::unique_ptr<Interior, InteriorDeleter>;
 
-TriangleMesh &      get_mesh(Interior &interior);
-const TriangleMesh &get_mesh(const Interior &interior);
+indexed_triangle_set &      get_mesh(Interior &interior);
+const indexed_triangle_set &get_mesh(const Interior &interior);
 
 struct DrainHole
 {
@@ -58,7 +56,7 @@ struct DrainHole
     bool get_intersections(const Vec3f& s, const Vec3f& dir,
                            std::array<std::pair<float, Vec3d>, 2>& out) const;
     
-    Contour3D to_mesh() const;
+    indexed_triangle_set to_mesh() const;
     
     template<class Archive> inline void serialize(Archive &ar)
     {
@@ -99,7 +97,13 @@ void cut_drainholes(std::vector<ExPolygons> & obj_slices,
                     const sla::DrainHoles &   holes,
                     std::function<void(void)> thr);
 
-}
+inline void swap_normals(indexed_triangle_set &its)
+{
+    for (auto &face : its.indices)
+        std::swap(face(0), face(2));
 }
 
+} // namespace sla
+} // namespace Slic3r
+
 #endif // HOLLOWINGFILTER_H
diff --git a/src/libslic3r/SLA/IndexedMesh.cpp b/src/libslic3r/SLA/IndexedMesh.cpp
index 2f47c63870..887ef1555d 100644
--- a/src/libslic3r/SLA/IndexedMesh.cpp
+++ b/src/libslic3r/SLA/IndexedMesh.cpp
@@ -10,59 +10,75 @@
 #include <libslic3r/SLA/Hollowing.hpp>
 #endif
 
-namespace Slic3r { namespace sla {
+namespace Slic3r {
+
+namespace sla {
 
 class IndexedMesh::AABBImpl {
 private:
     AABBTreeIndirect::Tree3f m_tree;
 
 public:
-    void init(const TriangleMesh& tm)
+    void init(const indexed_triangle_set &its)
     {
         m_tree = AABBTreeIndirect::build_aabb_tree_over_indexed_triangle_set(
-            tm.its.vertices, tm.its.indices);
+            its.vertices, its.indices);
     }
 
-    void intersect_ray(const TriangleMesh& tm,
-                       const Vec3d& s, const Vec3d& dir, igl::Hit& hit)
+    void intersect_ray(const indexed_triangle_set &its,
+                       const Vec3d &               s,
+                       const Vec3d &               dir,
+                       igl::Hit &                  hit)
     {
-        AABBTreeIndirect::intersect_ray_first_hit(tm.its.vertices,
-                                                  tm.its.indices,
-                                                  m_tree,
-                                                  s, dir, hit);
+        AABBTreeIndirect::intersect_ray_first_hit(its.vertices, its.indices,
+                                                  m_tree, s, dir, hit);
     }
 
-    void intersect_ray(const TriangleMesh& tm,
-                       const Vec3d& s, const Vec3d& dir, std::vector<igl::Hit>& hits)
+    void intersect_ray(const indexed_triangle_set &its,
+                       const Vec3d &               s,
+                       const Vec3d &               dir,
+                       std::vector<igl::Hit> &     hits)
     {
-        AABBTreeIndirect::intersect_ray_all_hits(tm.its.vertices,
-                                                 tm.its.indices,
-                                                 m_tree,
-                                                 s, dir, hits);
+        AABBTreeIndirect::intersect_ray_all_hits(its.vertices, its.indices,
+                                                 m_tree, s, dir, hits);
     }
 
-    double squared_distance(const TriangleMesh& tm,
-                            const Vec3d& point, int& i, Eigen::Matrix<double, 1, 3>& closest) {
+    double squared_distance(const indexed_triangle_set & its,
+                            const Vec3d &                point,
+                            int &                        i,
+                            Eigen::Matrix<double, 1, 3> &closest)
+    {
         size_t idx_unsigned = 0;
-        Vec3d closest_vec3d(closest);
-        double dist = AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
-            tm.its.vertices,
-            tm.its.indices,
-            m_tree, point, idx_unsigned, closest_vec3d);
-        i = int(idx_unsigned);
+        Vec3d  closest_vec3d(closest);
+        double dist =
+            AABBTreeIndirect::squared_distance_to_indexed_triangle_set(
+                its.vertices, its.indices, m_tree, point, idx_unsigned,
+                closest_vec3d);
+        i       = int(idx_unsigned);
         closest = closest_vec3d;
         return dist;
     }
 };
 
-IndexedMesh::IndexedMesh(const TriangleMesh& tmesh)
-    : m_aabb(new AABBImpl()), m_tm(&tmesh)
+template<class M> void IndexedMesh::init(const M &mesh)
 {
-    auto&& bb = tmesh.bounding_box();
+    BoundingBoxf3 bb = bounding_box(mesh);
     m_ground_level += bb.min(Z);
 
     // Build the AABB accelaration tree
-    m_aabb->init(tmesh);
+    m_aabb->init(*m_tm);
+}
+
+IndexedMesh::IndexedMesh(const indexed_triangle_set& tmesh)
+    : m_aabb(new AABBImpl()), m_tm(&tmesh)
+{
+    init(tmesh);
+}
+
+IndexedMesh::IndexedMesh(const TriangleMesh &mesh)
+    : m_aabb(new AABBImpl()), m_tm(&mesh.its)
+{
+    init(mesh);
 }
 
 IndexedMesh::~IndexedMesh() {}
@@ -87,34 +103,34 @@ IndexedMesh::IndexedMesh(IndexedMesh &&other) = default;
 
 const std::vector<Vec3f>& IndexedMesh::vertices() const
 {
-    return m_tm->its.vertices;
+    return m_tm->vertices;
 }
 
 
 
 const std::vector<Vec3i>& IndexedMesh::indices()  const
 {
-    return m_tm->its.indices;
+    return m_tm->indices;
 }
 
 
 
 const Vec3f& IndexedMesh::vertices(size_t idx) const
 {
-    return m_tm->its.vertices[idx];
+    return m_tm->vertices[idx];
 }
 
 
 
 const Vec3i& IndexedMesh::indices(size_t idx) const
 {
-    return m_tm->its.indices[idx];
+    return m_tm->indices[idx];
 }
 
 
-
 Vec3d IndexedMesh::normal_by_face_id(int face_id) const {
-    return m_tm->stl.facet_start[face_id].normal.cast<double>();
+
+    return its_unnormalized_normal(*m_tm, face_id).cast<double>().normalized();
 }
 
 
diff --git a/src/libslic3r/SLA/IndexedMesh.hpp b/src/libslic3r/SLA/IndexedMesh.hpp
index a72492b344..25ab75b836 100644
--- a/src/libslic3r/SLA/IndexedMesh.hpp
+++ b/src/libslic3r/SLA/IndexedMesh.hpp
@@ -15,6 +15,8 @@
   #include "libslic3r/SLA/Hollowing.hpp"
 #endif
 
+struct indexed_triangle_set;
+
 namespace Slic3r {
 
 class TriangleMesh;
@@ -29,7 +31,7 @@ using PointSet = Eigen::MatrixXd;
 class IndexedMesh {
     class AABBImpl;
     
-    const TriangleMesh* m_tm;
+    const indexed_triangle_set* m_tm;
     double m_ground_level = 0, m_gnd_offset = 0;
     
     std::unique_ptr<AABBImpl> m_aabb;
@@ -40,9 +42,12 @@ class IndexedMesh {
     std::vector<DrainHole> m_holes;
 #endif
 
+    template<class M> void init(const M &mesh);
+
 public:
     
-    explicit IndexedMesh(const TriangleMesh&);
+    explicit IndexedMesh(const indexed_triangle_set&);
+    explicit IndexedMesh(const TriangleMesh &mesh);
     
     IndexedMesh(const IndexedMesh& other);
     IndexedMesh& operator=(const IndexedMesh&);
@@ -130,7 +135,7 @@ public:
 
     Vec3d normal_by_face_id(int face_id) const;
 
-    const TriangleMesh * get_triangle_mesh() const { return m_tm; }
+    const indexed_triangle_set * get_triangle_mesh() const { return m_tm; }
 };
 
 // Calculate the normals for the selected points (from 'points' set) on the
diff --git a/src/libslic3r/SLA/Pad.cpp b/src/libslic3r/SLA/Pad.cpp
index bf2b4cf54a..8d995b59e3 100644
--- a/src/libslic3r/SLA/Pad.cpp
+++ b/src/libslic3r/SLA/Pad.cpp
@@ -1,7 +1,7 @@
 #include <libslic3r/SLA/Pad.hpp>
 #include <libslic3r/SLA/SpatIndex.hpp>
 #include <libslic3r/SLA/BoostAdapter.hpp>
-#include <libslic3r/SLA/Contour3D.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
 #include <libslic3r/TriangleMeshSlicer.hpp>
 
 #include "ConcaveHull.hpp"
@@ -29,25 +29,23 @@ namespace Slic3r { namespace sla {
 
 namespace {
 
-Contour3D walls(
+indexed_triangle_set walls(
     const Polygon &lower,
     const Polygon &upper,
     double         lower_z_mm,
     double         upper_z_mm)
 {
-    Wall w = triangulate_wall(lower, upper, lower_z_mm, upper_z_mm);
-
-    Contour3D ret;
-    ret.points = std::move(w.first);
-    ret.faces3 = std::move(w.second);
+    indexed_triangle_set w;
+    triangulate_wall(w.vertices, w.indices, lower, upper, lower_z_mm,
+                     upper_z_mm);
     
-    return ret;
+    return w;
 }
 
 // Same as walls() but with identical higher and lower polygons.
-Contour3D inline straight_walls(const Polygon &plate,
-                                double         lo_z,
-                                double         hi_z)
+inline indexed_triangle_set straight_walls(const Polygon &plate,
+                                           double         lo_z,
+                                           double         hi_z)
 {
     return walls(plate, plate, lo_z, hi_z);
 }
@@ -357,8 +355,10 @@ ExPolygon offset_contour_only(const ExPolygon &poly, coord_t delta, Args...args)
     return std::move(tmp2.front());
 }
 
-bool add_cavity(Contour3D &pad, ExPolygon &top_poly, const PadConfig3D &cfg,
-                ThrowOnCancel thr)
+bool add_cavity(indexed_triangle_set &pad,
+                ExPolygon &           top_poly,
+                const PadConfig3D &   cfg,
+                ThrowOnCancel         thr)
 {
     auto logerr = []{BOOST_LOG_TRIVIAL(error)<<"Could not create pad cavity";};
 
@@ -377,18 +377,18 @@ bool add_cavity(Contour3D &pad, ExPolygon &top_poly, const PadConfig3D &cfg,
     top_poly = pdiff.front();
 
     double z_min = -cfg.wing_height, z_max = 0;
-    pad.merge(walls(inner_base.contour, middle_base.contour, z_min, z_max));
+    its_merge(pad, walls(inner_base.contour, middle_base.contour, z_min, z_max));
     thr();
-    pad.merge(triangulate_expolygon_3d(inner_base, z_min, NORMALS_UP));
+    its_merge(pad, triangulate_expolygon_3d(inner_base, z_min, NORMALS_UP));
 
     return true;
 }
 
-Contour3D create_outer_pad_geometry(const ExPolygons & skeleton,
-                                    const PadConfig3D &cfg,
-                                    ThrowOnCancel      thr)
+indexed_triangle_set create_outer_pad_geometry(const ExPolygons & skeleton,
+                                               const PadConfig3D &cfg,
+                                               ThrowOnCancel      thr)
 {
-    Contour3D ret;
+    indexed_triangle_set ret;
 
     for (const ExPolygon &pad_part : skeleton) {
         ExPolygon top_poly{pad_part};
@@ -399,45 +399,45 @@ Contour3D create_outer_pad_geometry(const ExPolygons & skeleton,
         thr();
         
         double z_min = -cfg.height, z_max = 0;
-        ret.merge(walls(top_poly.contour, bottom_poly.contour, z_max, z_min));
+        its_merge(ret, walls(top_poly.contour, bottom_poly.contour, z_max, z_min));
 
         if (cfg.wing_height > 0. && add_cavity(ret, top_poly, cfg, thr))
             z_max = -cfg.wing_height;
 
         for (auto &h : bottom_poly.holes)
-            ret.merge(straight_walls(h, z_max, z_min));
+            its_merge(ret, straight_walls(h, z_max, z_min));
         
-        ret.merge(triangulate_expolygon_3d(bottom_poly, z_min, NORMALS_DOWN));
-        ret.merge(triangulate_expolygon_3d(top_poly, NORMALS_UP));
+        its_merge(ret, triangulate_expolygon_3d(bottom_poly, z_min, NORMALS_DOWN));
+        its_merge(ret, triangulate_expolygon_3d(top_poly, NORMALS_UP));
     }
 
     return ret;
 }
 
-Contour3D create_inner_pad_geometry(const ExPolygons & skeleton,
-                                    const PadConfig3D &cfg,
-                                    ThrowOnCancel      thr)
+indexed_triangle_set create_inner_pad_geometry(const ExPolygons & skeleton,
+                                               const PadConfig3D &cfg,
+                                               ThrowOnCancel      thr)
 {
-    Contour3D ret;
+    indexed_triangle_set ret;
 
     double z_max = 0., z_min = -cfg.height;
     for (const ExPolygon &pad_part : skeleton) {
         thr();
-        ret.merge(straight_walls(pad_part.contour, z_max, z_min));
+        its_merge(ret, straight_walls(pad_part.contour, z_max, z_min));
 
         for (auto &h : pad_part.holes)
-            ret.merge(straight_walls(h, z_max, z_min));
+            its_merge(ret, straight_walls(h, z_max, z_min));
     
-        ret.merge(triangulate_expolygon_3d(pad_part, z_min, NORMALS_DOWN));
-        ret.merge(triangulate_expolygon_3d(pad_part, z_max, NORMALS_UP));
+        its_merge(ret, triangulate_expolygon_3d(pad_part, z_min, NORMALS_DOWN));
+        its_merge(ret, triangulate_expolygon_3d(pad_part, z_max, NORMALS_UP));
     }
 
     return ret;
 }
 
-Contour3D create_pad_geometry(const PadSkeleton &skelet,
-                              const PadConfig &  cfg,
-                              ThrowOnCancel      thr)
+indexed_triangle_set create_pad_geometry(const PadSkeleton &skelet,
+                                         const PadConfig &  cfg,
+                                         ThrowOnCancel      thr)
 {
 #ifndef NDEBUG
     SVG svg("pad_skeleton.svg");
@@ -447,14 +447,16 @@ Contour3D create_pad_geometry(const PadSkeleton &skelet,
 #endif
 
     PadConfig3D cfg3d(cfg);
-    return create_outer_pad_geometry(skelet.outer, cfg3d, thr)
-        .merge(create_inner_pad_geometry(skelet.inner, cfg3d, thr));
+    auto pg = create_outer_pad_geometry(skelet.outer, cfg3d, thr);
+    its_merge(pg, create_inner_pad_geometry(skelet.inner, cfg3d, thr));
+
+    return pg;
 }
 
-Contour3D create_pad_geometry(const ExPolygons &supp_bp,
-                              const ExPolygons &model_bp,
-                              const PadConfig & cfg,
-                              ThrowOnCancel thr)
+indexed_triangle_set create_pad_geometry(const ExPolygons &supp_bp,
+                                         const ExPolygons &model_bp,
+                                         const PadConfig & cfg,
+                                         ThrowOnCancel     thr)
 {
     PadSkeleton skelet;
 
@@ -471,15 +473,14 @@ Contour3D create_pad_geometry(const ExPolygons &supp_bp,
 
 } // namespace
 
-void pad_blueprint(const TriangleMesh &      mesh,
-                   ExPolygons &              output,
-                   const std::vector<float> &heights,
-                   ThrowOnCancel             thrfn)
+void pad_blueprint(const indexed_triangle_set &mesh,
+                   ExPolygons &                output,
+                   const std::vector<float> &  heights,
+                   ThrowOnCancel               thrfn)
 {
     if (mesh.empty()) return;
 
-    assert(mesh.has_shared_vertices());
-    std::vector<ExPolygons> out = slice_mesh_ex(mesh.its, heights, thrfn);
+    std::vector<ExPolygons> out = slice_mesh_ex(mesh, heights, thrfn);
 
     size_t count = 0;
     for(auto& o : out) count += o.size();
@@ -500,26 +501,26 @@ void pad_blueprint(const TriangleMesh &      mesh,
     }
 }
 
-void pad_blueprint(const TriangleMesh &mesh,
-                   ExPolygons &        output,
-                   float               h,
-                   float               layerh,
-                   ThrowOnCancel       thrfn)
+void pad_blueprint(const indexed_triangle_set &mesh,
+                   ExPolygons &                output,
+                   float                       h,
+                   float                       layerh,
+                   ThrowOnCancel               thrfn)
 {
-    float gnd = float(mesh.bounding_box().min(Z));
+    float gnd = float(bounding_box(mesh).min(Z));
 
     std::vector<float> slicegrid = grid(gnd, gnd + h, layerh);
     pad_blueprint(mesh, output, slicegrid, thrfn);
 }
 
-void create_pad(const ExPolygons &sup_blueprint,
-                const ExPolygons &model_blueprint,
-                TriangleMesh &    out,
-                const PadConfig & cfg,
-                ThrowOnCancel thr)
+void create_pad(const ExPolygons &    sup_blueprint,
+                const ExPolygons &    model_blueprint,
+                indexed_triangle_set &out,
+                const PadConfig &     cfg,
+                ThrowOnCancel         thr)
 {
-    Contour3D t = create_pad_geometry(sup_blueprint, model_blueprint, cfg, thr);
-    out.merge(to_triangle_mesh(std::move(t)));
+    auto t = create_pad_geometry(sup_blueprint, model_blueprint, cfg, thr);
+    its_merge(out, t);
 }
 
 std::string PadConfig::validate() const
diff --git a/src/libslic3r/SLA/Pad.hpp b/src/libslic3r/SLA/Pad.hpp
index 61eec2dd1e..0b61495572 100644
--- a/src/libslic3r/SLA/Pad.hpp
+++ b/src/libslic3r/SLA/Pad.hpp
@@ -6,6 +6,8 @@
 #include <cmath>
 #include <string>
 
+struct indexed_triangle_set;
+
 namespace Slic3r {
 
 class ExPolygon;
@@ -13,25 +15,23 @@ class Polygon;
 using ExPolygons = std::vector<ExPolygon>;
 using Polygons = std::vector<Polygon>;
 
-class TriangleMesh;
-
 namespace sla {
 
 using ThrowOnCancel = std::function<void(void)>;
 
 /// Calculate the polygon representing the silhouette.
 void pad_blueprint(
-    const TriangleMesh &mesh,       // input mesh
+    const indexed_triangle_set &mesh,       // input mesh
     ExPolygons &        output,     // Output will be merged with
     const std::vector<float> &,     // Exact Z levels to sample
     ThrowOnCancel thrfn = [] {}); // Function that throws if cancel was requested
 
 void pad_blueprint(
-    const TriangleMesh &mesh,
-    ExPolygons &        output,
-    float               samplingheight = 0.1f,  // The height range to sample
-    float               layerheight    = 0.05f, // The sampling height
-    ThrowOnCancel       thrfn = [] {});
+    const indexed_triangle_set &mesh,
+    ExPolygons &                output,
+    float         samplingheight = 0.1f,  // The height range to sample
+    float         layerheight    = 0.05f, // The sampling height
+    ThrowOnCancel thrfn          = [] {});
 
 struct PadConfig {
     double wall_thickness_mm = 1.;
@@ -82,11 +82,12 @@ struct PadConfig {
     std::string validate() const;
 };
 
-void create_pad(const ExPolygons &support_contours,
-                const ExPolygons &model_contours,
-                TriangleMesh &    output_mesh,
-                const PadConfig & = PadConfig(),
-                ThrowOnCancel throw_on_cancel = []{});
+void create_pad(
+    const ExPolygons &    support_contours,
+    const ExPolygons &    model_contours,
+    indexed_triangle_set &output_mesh,
+    const PadConfig &             = PadConfig(),
+    ThrowOnCancel throw_on_cancel = [] {});
 
 } // namespace sla
 } // namespace Slic3r
diff --git a/src/libslic3r/SLA/SupportTree.cpp b/src/libslic3r/SLA/SupportTree.cpp
index 14a4dc3608..d9f8958a21 100644
--- a/src/libslic3r/SLA/SupportTree.cpp
+++ b/src/libslic3r/SLA/SupportTree.cpp
@@ -29,38 +29,36 @@
 namespace Slic3r {
 namespace sla {
 
-void SupportTree::retrieve_full_mesh(TriangleMesh &outmesh) const {
-    outmesh.merge(retrieve_mesh(MeshType::Support));
-    outmesh.merge(retrieve_mesh(MeshType::Pad));
+void SupportTree::retrieve_full_mesh(indexed_triangle_set &outmesh) const {
+    its_merge(outmesh, retrieve_mesh(MeshType::Support));
+    its_merge(outmesh, retrieve_mesh(MeshType::Pad));
 }
 
-std::vector<ExPolygons> SupportTree::slice(
-    const std::vector<float> &grid, float cr) const
+std::vector<ExPolygons> SupportTree::slice(const std::vector<float> &grid,
+                                           float                     cr) const
 {
-    const TriangleMesh &sup_mesh = retrieve_mesh(MeshType::Support);
-    const TriangleMesh &pad_mesh = retrieve_mesh(MeshType::Pad);
+    const indexed_triangle_set &sup_mesh = retrieve_mesh(MeshType::Support);
+    const indexed_triangle_set &pad_mesh = retrieve_mesh(MeshType::Pad);
 
     using Slices = std::vector<ExPolygons>;
     auto slices = reserve_vector<Slices>(2);
 
     if (!sup_mesh.empty()) {
         slices.emplace_back();
-        assert(sup_mesh.has_shared_vertices());
-        slices.back() = slice_mesh_ex(sup_mesh.its, grid, cr, ctl().cancelfn);
+        slices.back() = slice_mesh_ex(sup_mesh, grid, cr, ctl().cancelfn);
     }
 
     if (!pad_mesh.empty()) {
         slices.emplace_back();
 
-        auto bb = pad_mesh.bounding_box();
+        auto bb = bounding_box(pad_mesh);
         auto maxzit = std::upper_bound(grid.begin(), grid.end(), bb.max.z());
         
         auto cap = grid.end() - maxzit;
         auto padgrid = reserve_vector<float>(size_t(cap > 0 ? cap : 0));
         std::copy(grid.begin(), maxzit, std::back_inserter(padgrid));
 
-        assert(pad_mesh.has_shared_vertices());
-        slices.back() = slice_mesh_ex(pad_mesh.its, padgrid, cr, ctl().cancelfn);
+        slices.back() = slice_mesh_ex(pad_mesh, padgrid, cr, ctl().cancelfn);
     }
 
     size_t len = grid.size();
diff --git a/src/libslic3r/SLA/SupportTree.hpp b/src/libslic3r/SLA/SupportTree.hpp
index 4be90161d5..200b7cbde1 100644
--- a/src/libslic3r/SLA/SupportTree.hpp
+++ b/src/libslic3r/SLA/SupportTree.hpp
@@ -122,9 +122,9 @@ struct SupportableMesh
     IndexedMesh  emesh;
     SupportPoints pts;
     SupportTreeConfig cfg;
-    PadConfig     pad_cfg;
+//    PadConfig     pad_cfg;
 
-    explicit SupportableMesh(const TriangleMesh & trmsh,
+    explicit SupportableMesh(const indexed_triangle_set & trmsh,
                              const SupportPoints &sp,
                              const SupportTreeConfig &c)
         : emesh{trmsh}, pts{sp}, cfg{c}
@@ -149,22 +149,22 @@ public:
 
     virtual ~SupportTree() = default;
 
-    virtual const TriangleMesh &retrieve_mesh(MeshType meshtype) const = 0;
+    virtual const indexed_triangle_set &retrieve_mesh(MeshType meshtype) const = 0;
 
     /// Adding the "pad" under the supports.
     /// modelbase will be used according to the embed_object flag in PoolConfig.
     /// If set, the plate will be interpreted as the model's intrinsic pad. 
     /// Otherwise, the modelbase will be unified with the base plate calculated
     /// from the supports.
-    virtual const TriangleMesh &add_pad(const ExPolygons &modelbase,
-                                        const PadConfig & pcfg) = 0;
-    
+    virtual const indexed_triangle_set &add_pad(const ExPolygons &modelbase,
+                                                const PadConfig & pcfg) = 0;
+
     virtual void remove_pad() = 0;
     
     std::vector<ExPolygons> slice(const std::vector<float> &,
                                   float closing_radius) const;
     
-    void retrieve_full_mesh(TriangleMesh &outmesh) const;
+    void retrieve_full_mesh(indexed_triangle_set &outmesh) const;
     
     const JobController &ctl() const { return m_ctl; }
 };
diff --git a/src/libslic3r/SLA/SupportTreeBuilder.cpp b/src/libslic3r/SLA/SupportTreeBuilder.cpp
index daa01ef24d..86339d2acf 100644
--- a/src/libslic3r/SLA/SupportTreeBuilder.cpp
+++ b/src/libslic3r/SLA/SupportTreeBuilder.cpp
@@ -3,7 +3,7 @@
 #include <libslic3r/SLA/SupportTreeBuilder.hpp>
 #include <libslic3r/SLA/SupportTreeBuildsteps.hpp>
 #include <libslic3r/SLA/SupportTreeMesher.hpp>
-#include <libslic3r/SLA/Contour3D.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
 
 namespace Slic3r {
 namespace sla {
@@ -23,11 +23,11 @@ Head::Head(double       r_big_mm,
 {
 }
 
-Pad::Pad(const TriangleMesh &support_mesh,
-         const ExPolygons &  model_contours,
-         double              ground_level,
-         const PadConfig &   pcfg,
-         ThrowOnCancel       thr)
+Pad::Pad(const indexed_triangle_set &support_mesh,
+         const ExPolygons &          model_contours,
+         double                      ground_level,
+         const PadConfig &           pcfg,
+         ThrowOnCancel               thr)
     : cfg(pcfg)
     , zlevel(ground_level + pcfg.full_height() - pcfg.required_elevation())
 {
@@ -41,12 +41,14 @@ Pad::Pad(const TriangleMesh &support_mesh,
     pad_blueprint(support_mesh, sup_contours, grid(zstart, zend, 0.1f), thr);
     create_pad(sup_contours, model_contours, tmesh, pcfg);
     
-    tmesh.translate(0, 0, float(zlevel));
-    if (!tmesh.empty()) tmesh.require_shared_vertices();
+    Vec3f offs{.0f, .0f, float(zlevel)};
+    for (auto &p : tmesh.vertices) p += offs;
+
+    its_merge_vertices(tmesh);
 }
 
-const TriangleMesh &SupportTreeBuilder::add_pad(const ExPolygons &modelbase,
-                                                const PadConfig & cfg)
+const indexed_triangle_set &SupportTreeBuilder::add_pad(
+    const ExPolygons &modelbase, const PadConfig &cfg)
 {
     m_pad = Pad{merged_mesh(), modelbase, ground_level, cfg, ctl().cancelfn};
     return m_pad.tmesh;
@@ -120,74 +122,74 @@ void SupportTreeBuilder::add_pillar_base(long pid, double baseheight, double rad
     m_meshcache_valid = false;
 }
 
-const TriangleMesh &SupportTreeBuilder::merged_mesh(size_t steps) const
+const indexed_triangle_set &SupportTreeBuilder::merged_mesh(size_t steps) const
 {
     if (m_meshcache_valid) return m_meshcache;
     
-    Contour3D merged;
+    indexed_triangle_set merged;
     
     for (auto &head : m_heads) {
         if (ctl().stopcondition()) break;
-        if (head.is_valid()) merged.merge(get_mesh(head, steps));
+        if (head.is_valid()) its_merge(merged, get_mesh(head, steps));
     }
     
     for (auto &pill : m_pillars) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(pill, steps));
+        its_merge(merged, get_mesh(pill, steps));
     }
 
     for (auto &pedest : m_pedestals) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(pedest, steps));
+        its_merge(merged, get_mesh(pedest, steps));
     }
     
     for (auto &j : m_junctions) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(j, steps));
+        its_merge(merged, get_mesh(j, steps));
     }
 
     for (auto &bs : m_bridges) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(bs, steps));
+        its_merge(merged, get_mesh(bs, steps));
     }
     
     for (auto &bs : m_crossbridges) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(bs, steps));
+        its_merge(merged, get_mesh(bs, steps));
     }
 
     for (auto &bs : m_diffbridges) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(bs, steps));
+        its_merge(merged, get_mesh(bs, steps));
     }
 
     for (auto &anch : m_anchors) {
         if (ctl().stopcondition()) break;
-        merged.merge(get_mesh(anch, steps));
+        its_merge(merged, get_mesh(anch, steps));
     }
 
     if (ctl().stopcondition()) {
         // In case of failure we have to return an empty mesh
-        m_meshcache = TriangleMesh();
+        m_meshcache = {};
         return m_meshcache;
     }
     
-    m_meshcache = to_triangle_mesh(merged);
+    m_meshcache = merged;
     
     // The mesh will be passed by const-pointer to TriangleMeshSlicer,
     // which will need this.
-    if (!m_meshcache.empty()) m_meshcache.require_shared_vertices();
-    
-    BoundingBoxf3 &&bb = m_meshcache.bounding_box();
-    m_model_height       = bb.max(Z) - bb.min(Z);
+    its_merge_vertices(m_meshcache);
     
+    BoundingBoxf3 bb = bounding_box(m_meshcache);
+    m_model_height   = bb.max(Z) - bb.min(Z);
+
     m_meshcache_valid = true;
     return m_meshcache;
 }
 
 double SupportTreeBuilder::full_height() const
 {
-    if (merged_mesh().empty() && !pad().empty())
+    if (merged_mesh().indices.empty() && !pad().empty())
         return pad().cfg.full_height();
     
     double h = mesh_height();
@@ -195,7 +197,7 @@ double SupportTreeBuilder::full_height() const
     return h;
 }
 
-const TriangleMesh &SupportTreeBuilder::merge_and_cleanup()
+const indexed_triangle_set &SupportTreeBuilder::merge_and_cleanup()
 {
     // in case the mesh is not generated, it should be...
     auto &ret = merged_mesh(); 
@@ -210,7 +212,7 @@ const TriangleMesh &SupportTreeBuilder::merge_and_cleanup()
     return ret;
 }
 
-const TriangleMesh &SupportTreeBuilder::retrieve_mesh(MeshType meshtype) const
+const indexed_triangle_set &SupportTreeBuilder::retrieve_mesh(MeshType meshtype) const
 {
     switch(meshtype) {
     case MeshType::Support: return merged_mesh();
diff --git a/src/libslic3r/SLA/SupportTreeBuilder.hpp b/src/libslic3r/SLA/SupportTreeBuilder.hpp
index f29263ca3f..df5484e3dd 100644
--- a/src/libslic3r/SLA/SupportTreeBuilder.hpp
+++ b/src/libslic3r/SLA/SupportTreeBuilder.hpp
@@ -3,7 +3,8 @@
 
 #include <libslic3r/SLA/Concurrency.hpp>
 #include <libslic3r/SLA/SupportTree.hpp>
-#include <libslic3r/SLA/Contour3D.hpp>
+//#include <libslic3r/SLA/Contour3D.hpp>
+#include <libslic3r/TriangleMesh.hpp>
 #include <libslic3r/SLA/Pad.hpp>
 #include <libslic3r/MTUtils.hpp>
 
@@ -187,19 +188,19 @@ struct DiffBridge: public Bridge {
 
 // A wrapper struct around the pad
 struct Pad {
-    TriangleMesh tmesh;
+    indexed_triangle_set tmesh;
     PadConfig cfg;
     double zlevel = 0;
     
     Pad() = default;
-    
-    Pad(const TriangleMesh &support_mesh,
-        const ExPolygons &  model_contours,
-        double              ground_level,
-        const PadConfig &   pcfg,
-        ThrowOnCancel       thr);
-    
-    bool empty() const { return tmesh.facets_count() == 0; }
+
+    Pad(const indexed_triangle_set &support_mesh,
+        const ExPolygons &          model_contours,
+        double                      ground_level,
+        const PadConfig &           pcfg,
+        ThrowOnCancel               thr);
+
+    bool empty() const { return tmesh.indices.size() == 0; }
 };
 
 // This class will hold the support tree meshes with some additional
@@ -232,7 +233,7 @@ class SupportTreeBuilder: public SupportTree {
     
     using Mutex = ccr::SpinningMutex;
     
-    mutable TriangleMesh m_meshcache;
+    mutable indexed_triangle_set m_meshcache;
     mutable Mutex m_mutex;
     mutable bool m_meshcache_valid = false;
     mutable double m_model_height = 0; // the full height of the model
@@ -418,7 +419,7 @@ public:
     const Pad& pad() const { return m_pad; }
     
     // WITHOUT THE PAD!!!
-    const TriangleMesh &merged_mesh(size_t steps = 45) const;
+    const indexed_triangle_set &merged_mesh(size_t steps = 45) const;
     
     // WITH THE PAD
     double full_height() const;
@@ -431,16 +432,16 @@ public:
     }
     
     // Intended to be called after the generation is fully complete
-    const TriangleMesh & merge_and_cleanup();
+    const indexed_triangle_set & merge_and_cleanup();
     
     // Implement SupportTree interface:
 
-    const TriangleMesh &add_pad(const ExPolygons &modelbase,
-                                const PadConfig & pcfg) override;
-    
+    const indexed_triangle_set &add_pad(const ExPolygons &modelbase,
+                                        const PadConfig & pcfg) override;
+
     void remove_pad() override { m_pad = Pad(); }
-    
-    virtual const TriangleMesh &retrieve_mesh(
+
+    virtual const indexed_triangle_set &retrieve_mesh(
         MeshType meshtype = MeshType::Support) const override;
 };
 
diff --git a/src/libslic3r/SLA/SupportTreeMesher.cpp b/src/libslic3r/SLA/SupportTreeMesher.cpp
index 15491775b4..17f57f4578 100644
--- a/src/libslic3r/SLA/SupportTreeMesher.cpp
+++ b/src/libslic3r/SLA/SupportTreeMesher.cpp
@@ -2,22 +2,22 @@
 
 namespace Slic3r { namespace sla {
 
-Contour3D sphere(double rho, Portion portion, double fa) {
+indexed_triangle_set sphere(double rho, Portion portion, double fa) {
 
-    Contour3D ret;
+    indexed_triangle_set ret;
 
     // prohibit close to zero radius
     if(rho <= 1e-6 && rho >= -1e-6) return ret;
 
-    auto& vertices = ret.points;
-    auto& facets = ret.faces3;
+    auto& vertices = ret.vertices;
+    auto& facets = ret.indices;
 
     // Algorithm:
     // Add points one-by-one to the sphere grid and form facets using relative
     // coordinates. Sphere is composed effectively of a mesh of stacked circles.
 
     // adjust via rounding to get an even multiple for any provided angle.
-    double angle = (2*PI / floor(2*PI / fa));
+    double angle = (2 * PI / floor(2*PI / fa) );
 
     // Ring to be scaled to generate the steps of the sphere
     std::vector<double> ring;
@@ -32,8 +32,9 @@ Contour3D sphere(double rho, Portion portion, double fa) {
 
     // special case: first ring connects to 0,0,0
     // insert and form facets.
-    if(sbegin == 0)
-        vertices.emplace_back(Vec3d(0.0, 0.0, -rho + increment*sbegin*2.0*rho));
+    if (sbegin == 0)
+        vertices.emplace_back(
+            Vec3f(0.f, 0.f, float(-rho + increment * sbegin * 2. * rho)));
 
     auto id = coord_t(vertices.size());
     for (size_t i = 0; i < ring.size(); i++) {
@@ -42,7 +43,7 @@ Contour3D sphere(double rho, Portion portion, double fa) {
         // radius of the circle for this step.
         const double r = std::sqrt(std::abs(rho*rho - z*z));
         Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r);
-        vertices.emplace_back(Vec3d(b(0), b(1), z));
+        vertices.emplace_back(Vec3d(b(0), b(1), z).cast<float>());
 
         if (sbegin == 0)
             (i == 0) ? facets.emplace_back(coord_t(ring.size()), 0, 1) :
@@ -53,12 +54,12 @@ Contour3D sphere(double rho, Portion portion, double fa) {
     // General case: insert and form facets for each step,
     // joining it to the ring below it.
     for (size_t s = sbegin + 2; s < send - 1; s++) {
-        const double z = -rho + increment*double(s*2.0*rho);
+        const double z = -rho + increment * double(s * 2. * rho);
         const double r = std::sqrt(std::abs(rho*rho - z*z));
 
         for (size_t i = 0; i < ring.size(); i++) {
             Vec2d b = Eigen::Rotation2Dd(ring[i]) * Eigen::Vector2d(0, r);
-            vertices.emplace_back(Vec3d(b(0), b(1), z));
+            vertices.emplace_back(Vec3d(b(0), b(1), z).cast<float>());
             auto id_ringsize = coord_t(id - int(ring.size()));
             if (i == 0) {
                 // wrap around
@@ -75,7 +76,7 @@ Contour3D sphere(double rho, Portion portion, double fa) {
     // special case: last ring connects to 0,0,rho*2.0
     // only form facets.
     if(send >= size_t(2*PI / angle)) {
-        vertices.emplace_back(Vec3d(0.0, 0.0, -rho + increment*send*2.0*rho));
+        vertices.emplace_back(0.f, 0.f, float(-rho + increment*send*2.0*rho));
         for (size_t i = 0; i < ring.size(); i++) {
             auto id_ringsize = coord_t(id - int(ring.size()));
             if (i == 0) {
@@ -92,15 +93,15 @@ Contour3D sphere(double rho, Portion portion, double fa) {
     return ret;
 }
 
-Contour3D cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
+indexed_triangle_set cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
 {
     assert(ssteps > 0);
 
-    Contour3D ret;
+    indexed_triangle_set ret;
 
     auto steps = int(ssteps);
-    auto& points = ret.points;
-    auto& indices = ret.faces3;
+    auto& points = ret.vertices;
+    auto& indices = ret.indices;
     points.reserve(2*ssteps);
     double a = 2*PI/steps;
 
@@ -110,17 +111,17 @@ Contour3D cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
     // Upper circle points
     for(int i = 0; i < steps; ++i) {
         double phi = i*a;
-        double ex = endp(X) + r*std::cos(phi);
-        double ey = endp(Y) + r*std::sin(phi);
-        points.emplace_back(ex, ey, endp(Z));
+        auto ex = float(endp(X) + r*std::cos(phi));
+        auto ey = float(endp(Y) + r*std::sin(phi));
+        points.emplace_back(ex, ey, float(endp(Z)));
     }
 
     // Lower circle points
     for(int i = 0; i < steps; ++i) {
         double phi = i*a;
-        double x = jp(X) + r*std::cos(phi);
-        double y = jp(Y) + r*std::sin(phi);
-        points.emplace_back(x, y, jp(Z));
+        auto x = float(jp(X) + r*std::cos(phi));
+        auto y = float(jp(Y) + r*std::sin(phi));
+        points.emplace_back(x, y, float(jp(Z)));
     }
 
     // Now create long triangles connecting upper and lower circles
@@ -139,13 +140,13 @@ Contour3D cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
     // According to the slicing algorithms, we need to aid them with generating
     // a watertight body. So we create a triangle fan for the upper and lower
     // ending of the cylinder to close the geometry.
-    points.emplace_back(jp); int ci = int(points.size() - 1);
+    points.emplace_back(jp.cast<float>()); int ci = int(points.size() - 1);
     for(int i = 0; i < steps - 1; ++i)
         indices.emplace_back(i + offs + 1, i + offs, ci);
 
     indices.emplace_back(offs, steps + offs - 1, ci);
 
-    points.emplace_back(endp); ci = int(points.size() - 1);
+    points.emplace_back(endp.cast<float>()); ci = int(points.size() - 1);
     for(int i = 0; i < steps - 1; ++i)
         indices.emplace_back(ci, i, i + 1);
 
@@ -154,14 +155,17 @@ Contour3D cylinder(double r, double h, size_t ssteps, const Vec3d &sp)
     return ret;
 }
 
-Contour3D pinhead(double r_pin, double r_back, double length, size_t steps)
+indexed_triangle_set pinhead(double r_pin,
+                             double r_back,
+                             double length,
+                             size_t steps)
 {
     assert(steps > 0);
     assert(length >= 0.);
     assert(r_back > 0.);
     assert(r_pin > 0.);
 
-    Contour3D mesh;
+    indexed_triangle_set mesh;
 
     // We create two spheres which will be connected with a robe that fits
     // both circles perfectly.
@@ -187,66 +191,66 @@ Contour3D pinhead(double r_pin, double r_back, double length, size_t steps)
     auto &&s1 = sphere(r_back, make_portion(0, PI / 2 + phi), detail);
     auto &&s2 = sphere(r_pin, make_portion(PI / 2 + phi, PI), detail);
 
-    for (auto &p : s2.points) p.z() += h;
+    for (auto &p : s2.vertices) p.z() += h;
 
-    mesh.merge(s1);
-    mesh.merge(s2);
+    its_merge(mesh, s1);
+    its_merge(mesh, s2);
 
-    for (size_t idx1 = s1.points.size() - steps, idx2 = s1.points.size();
-         idx1 < s1.points.size() - 1; idx1++, idx2++) {
+    for (size_t idx1 = s1.vertices.size() - steps, idx2 = s1.vertices.size();
+         idx1 < s1.vertices.size() - 1; idx1++, idx2++) {
         coord_t i1s1 = coord_t(idx1), i1s2 = coord_t(idx2);
         coord_t i2s1 = i1s1 + 1, i2s2 = i1s2 + 1;
 
-        mesh.faces3.emplace_back(i1s1, i2s1, i2s2);
-        mesh.faces3.emplace_back(i1s1, i2s2, i1s2);
+        mesh.indices.emplace_back(i1s1, i2s1, i2s2);
+        mesh.indices.emplace_back(i1s1, i2s2, i1s2);
     }
 
-    auto i1s1 = coord_t(s1.points.size()) - coord_t(steps);
-    auto i2s1 = coord_t(s1.points.size()) - 1;
-    auto i1s2 = coord_t(s1.points.size());
-    auto i2s2 = coord_t(s1.points.size()) + coord_t(steps) - 1;
+    auto i1s1 = coord_t(s1.vertices.size()) - coord_t(steps);
+    auto i2s1 = coord_t(s1.vertices.size()) - 1;
+    auto i1s2 = coord_t(s1.vertices.size());
+    auto i2s2 = coord_t(s1.vertices.size()) + coord_t(steps) - 1;
 
-    mesh.faces3.emplace_back(i2s2, i2s1, i1s1);
-    mesh.faces3.emplace_back(i1s2, i2s2, i1s1);
+    mesh.indices.emplace_back(i2s2, i2s1, i1s1);
+    mesh.indices.emplace_back(i1s2, i2s2, i1s1);
 
     return mesh;
 }
 
-Contour3D halfcone(double       baseheight,
-                   double       r_bottom,
-                   double       r_top,
-                   const Vec3d &pos,
-                   size_t       steps)
+indexed_triangle_set halfcone(double       baseheight,
+                              double       r_bottom,
+                              double       r_top,
+                              const Vec3d &pos,
+                              size_t       steps)
 {
     assert(steps > 0);
 
     if (baseheight <= 0 || steps <= 0) return {};
 
-    Contour3D base;
+    indexed_triangle_set base;
 
     double a    = 2 * PI / steps;
     auto   last = int(steps - 1);
     Vec3d  ep{pos.x(), pos.y(), pos.z() + baseheight};
     for (size_t i = 0; i < steps; ++i) {
         double phi = i * a;
-        double x   = pos.x() + r_top * std::cos(phi);
-        double y   = pos.y() + r_top * std::sin(phi);
-        base.points.emplace_back(x, y, ep.z());
+        auto x   = float(pos.x() + r_top * std::cos(phi));
+        auto y   = float(pos.y() + r_top * std::sin(phi));
+        base.vertices.emplace_back(x, y, float(ep.z()));
     }
 
     for (size_t i = 0; i < steps; ++i) {
         double phi = i * a;
-        double x   = pos.x() + r_bottom * std::cos(phi);
-        double y   = pos.y() + r_bottom * std::sin(phi);
-        base.points.emplace_back(x, y, pos.z());
+        auto x   = float(pos.x() + r_bottom * std::cos(phi));
+        auto y   = float(pos.y() + r_bottom * std::sin(phi));
+        base.vertices.emplace_back(x, y, float(pos.z()));
     }
 
-    base.points.emplace_back(pos);
-    base.points.emplace_back(ep);
+    base.vertices.emplace_back(pos.cast<float>());
+    base.vertices.emplace_back(ep.cast<float>());
 
-    auto &indices = base.faces3;
-    auto  hcenter = int(base.points.size() - 1);
-    auto  lcenter = int(base.points.size() - 2);
+    auto &indices = base.indices;
+    auto  hcenter = int(base.vertices.size() - 1);
+    auto  lcenter = int(base.vertices.size() - 2);
     auto  offs    = int(steps);
     for (int i = 0; i < last; ++i) {
         indices.emplace_back(i, i + offs, offs + i + 1);
diff --git a/src/libslic3r/SLA/SupportTreeMesher.hpp b/src/libslic3r/SLA/SupportTreeMesher.hpp
index 63182745da..63ccc5fce3 100644
--- a/src/libslic3r/SLA/SupportTreeMesher.hpp
+++ b/src/libslic3r/SLA/SupportTreeMesher.hpp
@@ -4,7 +4,8 @@
 #include "libslic3r/Point.hpp"
 
 #include "libslic3r/SLA/SupportTreeBuilder.hpp"
-#include "libslic3r/SLA/Contour3D.hpp"
+#include "libslic3r/TriangleMesh.hpp"
+//#include "libslic3r/SLA/Contour3D.hpp"
 
 namespace Slic3r { namespace sla {
 
@@ -15,48 +16,53 @@ inline Portion make_portion(double a, double b)
     return std::make_tuple(a, b);
 }
 
-Contour3D sphere(double  rho,
-                 Portion portion = make_portion(0., 2. * PI),
-                 double  fa      = (2. * PI / 360.));
+indexed_triangle_set sphere(double  rho,
+                            Portion portion = make_portion(0., 2. * PI),
+                            double  fa      = (2. * PI / 360.));
 
 // Down facing cylinder in Z direction with arguments:
 // r: radius
 // h: Height
 // ssteps: how many edges will create the base circle
 // sp: starting point
-Contour3D cylinder(double       r,
-                   double       h,
-                   size_t       steps = 45,
-                   const Vec3d &sp    = Vec3d::Zero());
+indexed_triangle_set cylinder(double       r,
+                              double       h,
+                              size_t       steps = 45,
+                              const Vec3d &sp    = Vec3d::Zero());
 
-Contour3D pinhead(double r_pin, double r_back, double length, size_t steps = 45);
+indexed_triangle_set pinhead(double r_pin,
+                             double r_back,
+                             double length,
+                             size_t steps = 45);
 
-Contour3D halfcone(double       baseheight,
-                   double       r_bottom,
-                   double       r_top,
-                   const Vec3d &pt    = Vec3d::Zero(),
-                   size_t       steps = 45);
+indexed_triangle_set halfcone(double       baseheight,
+                              double       r_bottom,
+                              double       r_top,
+                              const Vec3d &pt    = Vec3d::Zero(),
+                              size_t       steps = 45);
 
-inline Contour3D get_mesh(const Head &h, size_t steps)
+inline indexed_triangle_set get_mesh(const Head &h, size_t steps)
 {
-    Contour3D mesh = pinhead(h.r_pin_mm, h.r_back_mm, h.width_mm, steps);
+    indexed_triangle_set mesh = pinhead(h.r_pin_mm, h.r_back_mm, h.width_mm, steps);
 
-    for(auto& p : mesh.points) p.z() -= (h.fullwidth() - h.r_back_mm);
+    for (auto& p : mesh.vertices) p.z() -= (h.fullwidth() - h.r_back_mm);
 
-    using Quaternion = Eigen::Quaternion<double>;
+    using Quaternion = Eigen::Quaternion<float>;
 
     // We rotate the head to the specified direction. The head's pointing
     // side is facing upwards so this means that it would hold a support
     // point with a normal pointing straight down. This is the reason of
     // the -1 z coordinate
-    auto quatern = Quaternion::FromTwoVectors(Vec3d{0, 0, -1}, h.dir);
+    auto quatern = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, -1.f},
+                                              h.dir.cast<float>());
 
-    for(auto& p : mesh.points) p = quatern * p + h.pos;
+    Vec3f pos = h.pos.cast<float>();
+    for (auto& p : mesh.vertices) p = quatern * p + pos;
 
     return mesh;
 }
 
-inline Contour3D get_mesh(const Pillar &p, size_t steps)
+inline indexed_triangle_set get_mesh(const Pillar &p, size_t steps)
 {
     if(p.height > EPSILON) { // Endpoint is below the starting point
         // We just create a bridge geometry with the pillar parameters and
@@ -67,47 +73,53 @@ inline Contour3D get_mesh(const Pillar &p, size_t steps)
     return {};
 }
 
-inline Contour3D get_mesh(const Pedestal &p, size_t steps)
+inline indexed_triangle_set get_mesh(const Pedestal &p, size_t steps)
 {
     return halfcone(p.height, p.r_bottom, p.r_top, p.pos, steps);
 }
 
-inline Contour3D get_mesh(const Junction &j, size_t steps)
+inline indexed_triangle_set get_mesh(const Junction &j, size_t steps)
 {
-    Contour3D mesh = sphere(j.r, make_portion(0, PI), 2 *PI / steps);
-    for(auto& p : mesh.points) p += j.pos;
+    indexed_triangle_set mesh = sphere(j.r, make_portion(0, PI), 2 *PI / steps);
+    Vec3f pos = j.pos.cast<float>();
+    for(auto& p : mesh.vertices) p += pos;
     return mesh;
 }
 
-inline Contour3D get_mesh(const Bridge &br, size_t steps)
+inline indexed_triangle_set get_mesh(const Bridge &br, size_t steps)
 {
-    using Quaternion = Eigen::Quaternion<double>;
+    using Quaternion = Eigen::Quaternion<float>;
     Vec3d v = (br.endp - br.startp);
     Vec3d dir = v.normalized();
     double d = v.norm();
 
-    Contour3D mesh = cylinder(br.r, d, steps);
+    indexed_triangle_set mesh = cylinder(br.r, d, steps);
 
-    auto quater = Quaternion::FromTwoVectors(Vec3d{0,0,1}, dir);
-    for(auto& p : mesh.points) p = quater * p + br.startp;
+    auto quater = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f},
+                                             dir.cast<float>());
+
+    Vec3f startp = br.startp.cast<float>();
+    for(auto& p : mesh.vertices) p = quater * p + startp;
 
     return mesh;
 }
 
-inline Contour3D get_mesh(const DiffBridge &br, size_t steps)
+inline indexed_triangle_set get_mesh(const DiffBridge &br, size_t steps)
 {
     double h = br.get_length();
-    Contour3D mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps);
+    indexed_triangle_set mesh = halfcone(h, br.r, br.end_r, Vec3d::Zero(), steps);
 
-    using Quaternion = Eigen::Quaternion<double>;
+    using Quaternion = Eigen::Quaternion<float>;
 
     // We rotate the head to the specified direction. The head's pointing
     // side is facing upwards so this means that it would hold a support
     // point with a normal pointing straight down. This is the reason of
     // the -1 z coordinate
-    auto quatern = Quaternion::FromTwoVectors(Vec3d{0, 0, 1}, br.get_dir());
+    auto quatern = Quaternion::FromTwoVectors(Vec3f{0.f, 0.f, 1.f},
+                                              br.get_dir().cast<float>());
 
-    for(auto& p : mesh.points) p = quatern * p + br.startp;
+    Vec3f startp = br.startp.cast<float>();
+    for(auto& p : mesh.vertices) p = quatern * p + startp;
 
     return mesh;
 }
diff --git a/src/libslic3r/SLAPrint.cpp b/src/libslic3r/SLAPrint.cpp
index c0c65d39bf..27423793ee 100644
--- a/src/libslic3r/SLAPrint.cpp
+++ b/src/libslic3r/SLAPrint.cpp
@@ -109,7 +109,7 @@ sla::PadConfig make_pad_cfg(const SLAPrintObjectConfig& c)
     return pcfg;
 }
 
-bool validate_pad(const TriangleMesh &pad, const sla::PadConfig &pcfg) 
+bool validate_pad(const indexed_triangle_set &pad, const sla::PadConfig &pcfg)
 {
     // An empty pad can only be created if embed_object mode is enabled
     // and the pad is not forced everywhere
@@ -1067,6 +1067,7 @@ Vec3d SLAPrint::relative_correction() const
 namespace { // dummy empty static containers for return values in some methods
 const std::vector<ExPolygons> EMPTY_SLICES;
 const TriangleMesh EMPTY_MESH;
+const indexed_triangle_set EMPTY_TRIANGLE_SET;
 const ExPolygons EMPTY_SLICE;
 const std::vector<sla::SupportPoint> EMPTY_SUPPORT_POINTS;
 }
@@ -1129,31 +1130,27 @@ TriangleMesh SLAPrintObject::get_mesh(SLAPrintObjectStep step) const
 
 const TriangleMesh& SLAPrintObject::support_mesh() const
 {
-    sla::SupportTree::UPtr &stree = m_supportdata->support_tree_ptr;
-    
-    if(m_config.supports_enable.getBool() && m_supportdata && stree)
-        return stree->retrieve_mesh(sla::MeshType::Support);
+    if(m_config.supports_enable.getBool() && m_supportdata)
+        return m_supportdata->tree_mesh;
     
     return EMPTY_MESH;
 }
 
 const TriangleMesh& SLAPrintObject::pad_mesh() const
 {
-    sla::SupportTree::UPtr &stree = m_supportdata->support_tree_ptr;
-    
-    if(m_config.pad_enable.getBool() && m_supportdata && stree)
-        return stree->retrieve_mesh(sla::MeshType::Pad);
+    if(m_config.pad_enable.getBool() && m_supportdata)
+        return m_supportdata->pad_mesh;
 
     return EMPTY_MESH;
 }
 
-const TriangleMesh &SLAPrintObject::hollowed_interior_mesh() const
+const indexed_triangle_set &SLAPrintObject::hollowed_interior_mesh() const
 {
     if (m_hollowing_data && m_hollowing_data->interior &&
         m_config.hollowing_enable.getBool())
         return sla::get_mesh(*m_hollowing_data->interior);
     
-    return EMPTY_MESH;
+    return EMPTY_TRIANGLE_SET;
 }
 
 const TriangleMesh &SLAPrintObject::transformed_mesh() const {
diff --git a/src/libslic3r/SLAPrint.hpp b/src/libslic3r/SLAPrint.hpp
index f5f422b3d8..e11926c7ea 100644
--- a/src/libslic3r/SLAPrint.hpp
+++ b/src/libslic3r/SLAPrint.hpp
@@ -79,8 +79,8 @@ public:
     const TriangleMesh&     pad_mesh() const;
     
     // Ready after this->is_step_done(slaposDrillHoles) is true
-    const TriangleMesh&     hollowed_interior_mesh() const;
-    
+    const indexed_triangle_set &hollowed_interior_mesh() const;
+
     // Get the mesh that is going to be printed with all the modifications
     // like hollowing and drilled holes.
     const TriangleMesh & get_mesh_to_print() const {
@@ -313,16 +313,27 @@ private:
     public:
         sla::SupportTree::UPtr  support_tree_ptr; // the supports
         std::vector<ExPolygons> support_slices;   // sliced supports
+        TriangleMesh tree_mesh, pad_mesh, full_mesh;
         
         inline SupportData(const TriangleMesh &t)
-            : sla::SupportableMesh{t, {}, {}}
+            : sla::SupportableMesh{t.its, {}, {}}
         {}
         
         sla::SupportTree::UPtr &create_support_tree(const sla::JobController &ctl)
         {
             support_tree_ptr = sla::SupportTree::create(*this, ctl);
+            tree_mesh = TriangleMesh{support_tree_ptr->retrieve_mesh(sla::MeshType::Support)};
             return support_tree_ptr;
         }
+
+        void create_pad(const ExPolygons &blueprint, const sla::PadConfig &pcfg)
+        {
+            if (!support_tree_ptr)
+                return;
+
+            support_tree_ptr->add_pad(blueprint, pcfg);
+            pad_mesh = TriangleMesh{support_tree_ptr->retrieve_mesh(sla::MeshType::Pad)};
+        }
     };
     
     std::unique_ptr<SupportData> m_supportdata;
@@ -569,7 +580,7 @@ sla::PadConfig::EmbedObject builtin_pad_cfg(const SLAPrintObjectConfig& c);
 
 sla::PadConfig make_pad_cfg(const SLAPrintObjectConfig& c);
 
-bool validate_pad(const TriangleMesh &pad, const sla::PadConfig &pcfg);
+bool validate_pad(const indexed_triangle_set &pad, const sla::PadConfig &pcfg);
 
 
 } // namespace Slic3r
diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp
index 94b4c57a42..d6620b453f 100644
--- a/src/libslic3r/SLAPrintSteps.cpp
+++ b/src/libslic3r/SLAPrintSteps.cpp
@@ -194,18 +194,18 @@ static std::vector<bool> create_exclude_mask(
         const sla::Interior &interior,
         const std::vector<sla::DrainHole> &holes)
 {
-    FaceHash interior_hash{sla::get_mesh(interior).its};
+    FaceHash interior_hash{sla::get_mesh(interior)};
 
     std::vector<bool> exclude_mask(its.indices.size(), false);
 
-    std::vector< std::vector<size_t> > neighbor_index =
-            create_vertex_faces_index(its);
+    VertexFaceIndex neighbor_index{its};
 
     auto exclude_neighbors = [&neighbor_index, &exclude_mask](const Vec3i &face)
     {
         for (int i = 0; i < 3; ++i) {
-            const std::vector<size_t> &neighbors = neighbor_index[face(i)];
-            for (size_t fi_n : neighbors) exclude_mask[fi_n] = true;
+            const auto &neighbors_range = neighbor_index[face(i)];
+            for (size_t fi_n : neighbors_range)
+                exclude_mask[fi_n] = true;
         }
     };
 
@@ -360,7 +360,7 @@ void SLAPrint::Steps::drill_holes(SLAPrintObject &po)
         holept.normal += Vec3f{dist(m_rng), dist(m_rng), dist(m_rng)};
         holept.normal.normalize();
         holept.pos += Vec3f{dist(m_rng), dist(m_rng), dist(m_rng)};
-        TriangleMesh m = sla::to_triangle_mesh(holept.to_mesh());
+        TriangleMesh m{holept.to_mesh()};
         m.require_shared_vertices();
 
         part_to_drill.indices.clear();
@@ -489,11 +489,11 @@ void SLAPrint::Steps::slice_model(SLAPrintObject &po)
                                   nullptr;
 
     if (interior && ! sla::get_mesh(*interior).empty()) {
-        TriangleMesh interiormesh = sla::get_mesh(*interior);
-        interiormesh.repaired = false;
-        interiormesh.repair(true);
+        indexed_triangle_set interiormesh = sla::get_mesh(*interior);
+        sla::swap_normals(interiormesh);
         params.mode = MeshSlicingParams::SlicingMode::Regular;
-        std::vector<ExPolygons> interior_slices = slice_mesh_ex(interiormesh.its, slice_grid, params, thr);
+
+        std::vector<ExPolygons> interior_slices = slice_mesh_ex(interiormesh, slice_grid, params, thr);
 
         sla::ccr::for_each(size_t(0), interior_slices.size(),
                            [&po, &interior_slices] (size_t i) {
@@ -667,15 +667,14 @@ void SLAPrint::Steps::generate_pad(SLAPrintObject &po) {
             // we sometimes call it "builtin pad" is enabled so we will
             // get a sample from the bottom of the mesh and use it for pad
             // creation.
-            sla::pad_blueprint(trmesh, bp, float(pad_h),
+            sla::pad_blueprint(trmesh.its, bp, float(pad_h),
                                float(po.m_config.layer_height.getFloat()),
                                [this](){ throw_if_canceled(); });
         }
 
-        po.m_supportdata->support_tree_ptr->add_pad(bp, pcfg);
-        auto &pad_mesh = po.m_supportdata->support_tree_ptr->retrieve_mesh(sla::MeshType::Pad);
+        po.m_supportdata->create_pad(bp, pcfg);
 
-        if (!validate_pad(pad_mesh, pcfg))
+        if (!validate_pad(po.m_supportdata->support_tree_ptr->retrieve_mesh(sla::MeshType::Pad), pcfg))
             throw Slic3r::SlicingError(
                     L("No pad can be generated for this model with the "
                       "current configuration"));
diff --git a/src/libslic3r/SlicesToTriangleMesh.cpp b/src/libslic3r/SlicesToTriangleMesh.cpp
index d6a5469619..7a2975d12d 100644
--- a/src/libslic3r/SlicesToTriangleMesh.cpp
+++ b/src/libslic3r/SlicesToTriangleMesh.cpp
@@ -1,8 +1,8 @@
 
 #include "SlicesToTriangleMesh.hpp"
 
-#include "libslic3r/MTUtils.hpp"
-#include "libslic3r/SLA/Contour3D.hpp"
+//#include "libslic3r/MTUtils.hpp"
+#include "libslic3r/Execution/ExecutionTBB.hpp"
 #include "libslic3r/ClipperUtils.hpp"
 #include "libslic3r/Tesselate.hpp"
 
@@ -11,71 +11,73 @@
 
 namespace Slic3r {
 
-inline sla::Contour3D wall_strip(const Polygon &poly,
-                                 double         lower_z_mm,
-                                 double         upper_z_mm)
+inline indexed_triangle_set wall_strip(const Polygon &poly,
+                                       double         lower_z_mm,
+                                       double         upper_z_mm)
 {
-    sla::Contour3D ret;
+    indexed_triangle_set ret;
     
-    size_t startidx = ret.points.size();
+    size_t startidx = ret.vertices.size();
     size_t offs     = poly.points.size();
     
-    ret.points.reserve(ret.points.size() + 2 *offs);
+    ret.vertices.reserve(ret.vertices.size() + 2 *offs);
     
     for (const Point &p : poly.points)
-        ret.points.emplace_back(to_3d(unscaled(p), lower_z_mm));
+        ret.vertices.emplace_back(to_3d(unscaled<float>(p), float(lower_z_mm)));
     
     for (const Point &p : poly.points)
-        ret.points.emplace_back(to_3d(unscaled(p), upper_z_mm));
+        ret.vertices.emplace_back(to_3d(unscaled<float>(p), float(upper_z_mm)));
     
     for (size_t i = startidx + 1; i < startidx + offs; ++i) {
-        ret.faces3.emplace_back(i - 1, i, i + offs - 1);
-        ret.faces3.emplace_back(i, i + offs, i + offs - 1);
+        ret.indices.emplace_back(i - 1, i, i + offs - 1);
+        ret.indices.emplace_back(i, i + offs, i + offs - 1);
     }
     
-    ret.faces3.emplace_back(startidx + offs - 1, startidx, startidx + 2 * offs - 1);
-    ret.faces3.emplace_back(startidx, startidx + offs, startidx + 2 * offs - 1);
+    ret.indices.emplace_back(startidx + offs - 1, startidx, startidx + 2 * offs - 1);
+    ret.indices.emplace_back(startidx, startidx + offs, startidx + 2 * offs - 1);
     
     return ret;
 }
 
 // Same as walls() but with identical higher and lower polygons.
-sla::Contour3D inline straight_walls(const Polygon &plate,
+indexed_triangle_set inline straight_walls(const Polygon &plate,
                                      double         lo_z,
                                      double         hi_z)
 {
     return wall_strip(plate, lo_z, hi_z);
 }
 
-sla::Contour3D inline straight_walls(const ExPolygon &plate,
+indexed_triangle_set inline straight_walls(const ExPolygon &plate,
                                      double           lo_z,
                                      double           hi_z)
 {
-    sla::Contour3D ret;
-    ret.merge(straight_walls(plate.contour, lo_z, hi_z));
-    for (auto &h : plate.holes) ret.merge(straight_walls(h, lo_z, hi_z));
+    indexed_triangle_set ret = straight_walls(plate.contour, lo_z, hi_z);
+    for (auto &h : plate.holes)
+        its_merge(ret, straight_walls(h, lo_z, hi_z));
+
     return ret;
 }
 
-sla::Contour3D inline straight_walls(const ExPolygons &slice,
+indexed_triangle_set inline straight_walls(const ExPolygons &slice,
                                      double            lo_z,
                                      double            hi_z)
-{   
-    sla::Contour3D ret;
+{
+    indexed_triangle_set ret;
     for (const ExPolygon &poly : slice)
-        ret.merge(straight_walls(poly, lo_z, hi_z));
-    
+        its_merge(ret, straight_walls(poly, lo_z, hi_z));
+
     return ret;
 }
 
-sla::Contour3D slices_to_triangle_mesh(const std::vector<ExPolygons> &slices,
-                                       double zmin,
-                                       const std::vector<float> &    grid)
+indexed_triangle_set slices_to_mesh(
+    const std::vector<ExPolygons> &slices,
+    double                         zmin,
+    const std::vector<float> &     grid)
 {
     assert(slices.size() == grid.size());
 
-    using Layers = std::vector<sla::Contour3D>;
-    std::vector<sla::Contour3D> layers(slices.size());
+    using Layers = std::vector<indexed_triangle_set>;
+    Layers layers(slices.size());
     size_t len = slices.size() - 1;
 
     tbb::parallel_for(size_t(0), len, [&slices, &layers, &grid](size_t i) {
@@ -84,45 +86,52 @@ sla::Contour3D slices_to_triangle_mesh(const std::vector<ExPolygons> &slices,
 
         ExPolygons dff1 = diff_ex(lower, upper);
         ExPolygons dff2 = diff_ex(upper, lower);
-        layers[i].merge(triangulate_expolygons_3d(dff1, grid[i], NORMALS_UP));
-        layers[i].merge(triangulate_expolygons_3d(dff2, grid[i], NORMALS_DOWN));
-        layers[i].merge(straight_walls(upper, grid[i], grid[i + 1]));
+        its_merge(layers[i], triangulate_expolygons_3d(dff1, grid[i], NORMALS_UP));
+        its_merge(layers[i], triangulate_expolygons_3d(dff2, grid[i], NORMALS_DOWN));
+        its_merge(layers[i], straight_walls(upper, grid[i], grid[i + 1]));
         
     });
+
+    auto merge_fn = []( const indexed_triangle_set &a, const indexed_triangle_set &b ) {
+        indexed_triangle_set res{a}; its_merge(res, b); return res;
+    };
+
+    auto ret = execution::reduce(ex_tbb, layers.begin(), layers.end(),
+                                 indexed_triangle_set{}, merge_fn);
+
+    //    sla::Contour3D ret = tbb::parallel_reduce(
+    //        tbb::blocked_range(layers.begin(), layers.end()),
+    //        sla::Contour3D{},
+    //        [](const tbb::blocked_range<Layers::iterator>& r, sla::Contour3D
+    //        init) {
+    //            for(auto it = r.begin(); it != r.end(); ++it )
+    //            init.merge(*it); return init;
+    //        },
+    //        []( const sla::Contour3D &a, const sla::Contour3D &b ) {
+    //            sla::Contour3D res{a}; res.merge(b); return res;
+    //        });
     
-    sla::Contour3D ret = tbb::parallel_reduce(
-        tbb::blocked_range(layers.begin(), layers.end()),
-        sla::Contour3D{},
-        [](const tbb::blocked_range<Layers::iterator>& r, sla::Contour3D init) {
-            for(auto it = r.begin(); it != r.end(); ++it ) init.merge(*it);
-            return init;
-        },
-        []( const sla::Contour3D &a, const sla::Contour3D &b ) {
-            sla::Contour3D res{a}; res.merge(b); return res;
-        });
-    
-    ret.merge(triangulate_expolygons_3d(slices.front(), zmin, NORMALS_DOWN));
-    ret.merge(straight_walls(slices.front(), zmin, grid.front()));
-    ret.merge(triangulate_expolygons_3d(slices.back(), grid.back(), NORMALS_UP));
+    its_merge(ret, triangulate_expolygons_3d(slices.front(), zmin, NORMALS_DOWN));
+    its_merge(ret, straight_walls(slices.front(), zmin, grid.front()));
+    its_merge(ret, triangulate_expolygons_3d(slices.back(), grid.back(), NORMALS_UP));
         
     return ret;
 }
 
-void slices_to_triangle_mesh(TriangleMesh &                 mesh,
+void slices_to_mesh(indexed_triangle_set &         mesh,
                              const std::vector<ExPolygons> &slices,
                              double                         zmin,
                              double                         lh,
                              double                         ilh)
 {
-    std::vector<sla::Contour3D> wall_meshes(slices.size());
+    std::vector<indexed_triangle_set> wall_meshes(slices.size());
     std::vector<float> grid(slices.size(), zmin + ilh);
     
-    for (size_t i = 1; i < grid.size(); ++i) grid[i] = grid[i - 1] + lh;
+    for (size_t i = 1; i < grid.size(); ++i)
+        grid[i] = grid[i - 1] + lh;
     
-    sla::Contour3D cntr = slices_to_triangle_mesh(slices, zmin, grid);
-    mesh.merge(sla::to_triangle_mesh(cntr));
-    mesh.repaired = true;
-    mesh.require_shared_vertices();
+    indexed_triangle_set cntr = slices_to_mesh(slices, zmin, grid);
+    its_merge(mesh, cntr);
 }
 
 } // namespace Slic3r
diff --git a/src/libslic3r/SlicesToTriangleMesh.hpp b/src/libslic3r/SlicesToTriangleMesh.hpp
index 133312d569..2fd1778851 100644
--- a/src/libslic3r/SlicesToTriangleMesh.hpp
+++ b/src/libslic3r/SlicesToTriangleMesh.hpp
@@ -6,16 +6,18 @@
 
 namespace Slic3r {
 
-void slices_to_triangle_mesh(TriangleMesh &                 mesh,
+void slices_to_mesh(indexed_triangle_set &         mesh,
                              const std::vector<ExPolygons> &slices,
                              double                         zmin,
                              double                         lh,
                              double                         ilh);
 
-inline TriangleMesh slices_to_triangle_mesh(
+inline indexed_triangle_set slices_to_mesh(
     const std::vector<ExPolygons> &slices, double zmin, double lh, double ilh)
 {
-    TriangleMesh out; slices_to_triangle_mesh(out, slices, zmin, lh, ilh);
+    indexed_triangle_set out;
+    slices_to_mesh(out, slices, zmin, lh, ilh);
+
     return out;
 }
 
diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp
index a8e952a2b2..93276d6006 100644
--- a/src/libslic3r/TriangleMesh.cpp
+++ b/src/libslic3r/TriangleMesh.cpp
@@ -1,8 +1,12 @@
 #include "Exception.hpp"
 #include "TriangleMesh.hpp"
 #include "TriangleMeshSlicer.hpp"
+#include "MeshSplitImpl.hpp"
 #include "ClipperUtils.hpp"
 #include "Geometry.hpp"
+#include "Point.hpp"
+#include "Execution/ExecutionTBB.hpp"
+#include "Execution/ExecutionSeq.hpp"
 
 #include <libqhullcpp/Qhull.h>
 #include <libqhullcpp/QhullFacetList.h>
@@ -447,27 +451,40 @@ std::deque<uint32_t> TriangleMesh::find_unvisited_neighbors(std::vector<unsigned
  */
 TriangleMeshPtrs TriangleMesh::split() const
 {
-    // Loop while we have remaining facets.
-    std::vector<unsigned char> facet_visited;
+    struct MeshAdder {
+        TriangleMeshPtrs &meshes;
+        MeshAdder(TriangleMeshPtrs &ptrs): meshes{ptrs} {}
+        void operator=(const indexed_triangle_set &its)
+        {
+            meshes.emplace_back(new TriangleMesh(its));
+        }
+    };
+
     TriangleMeshPtrs meshes;
-    for (;;) {
-        std::deque<uint32_t> facets = find_unvisited_neighbors(facet_visited);
-        if (facets.empty())
-            break;
+    if (has_shared_vertices()) {
+        its_split(its, MeshAdder{meshes});
+    } else {
+        // Loop while we have remaining facets.
+        std::vector<unsigned char> facet_visited;
+        for (;;) {
+            std::deque<uint32_t> facets = find_unvisited_neighbors(facet_visited);
+            if (facets.empty())
+                break;
 
-        // Create a new mesh for the part that was just split off.
-        TriangleMesh* mesh = new TriangleMesh;
-        meshes.emplace_back(mesh);
-        mesh->stl.stats.type = inmemory;
-        mesh->stl.stats.number_of_facets = (uint32_t)facets.size();
-        mesh->stl.stats.original_num_facets = mesh->stl.stats.number_of_facets;
-        stl_allocate(&mesh->stl);
+            // Create a new mesh for the part that was just split off.
+            TriangleMesh* mesh = new TriangleMesh;
+            meshes.emplace_back(mesh);
+            mesh->stl.stats.type = inmemory;
+            mesh->stl.stats.number_of_facets = (uint32_t)facets.size();
+            mesh->stl.stats.original_num_facets = mesh->stl.stats.number_of_facets;
+            stl_allocate(&mesh->stl);
 
-        // Assign the facets to the new mesh.
-        bool first = true;
-        for (auto facet = facets.begin(); facet != facets.end(); ++ facet) {
-            mesh->stl.facet_start[facet - facets.begin()] = this->stl.facet_start[*facet];
-            stl_facet_stats(&mesh->stl, this->stl.facet_start[*facet], first);
+            // Assign the facets to the new mesh.
+            bool first = true;
+            for (auto facet = facets.begin(); facet != facets.end(); ++ facet) {
+                mesh->stl.facet_start[facet - facets.begin()] = this->stl.facet_start[*facet];
+                stl_facet_stats(&mesh->stl, this->stl.facet_start[*facet], first);
+            }
         }
     }
 
@@ -667,69 +684,24 @@ void TriangleMesh::restore_optional()
     }
 }
 
-std::vector<std::vector<size_t>> create_vertex_faces_index(const indexed_triangle_set &its)
-{
-    std::vector<std::vector<size_t>> index;
+// Create a mapping from triangle edge into face.
+struct EdgeToFace {
+    // Index of the 1st vertex of the triangle edge. vertex_low <= vertex_high.
+    int  vertex_low;
+    // Index of the 2nd vertex of the triangle edge.
+    int  vertex_high;
+    // Index of a triangular face.
+    int  face;
+    // Index of edge in the face, starting with 1. Negative indices if the edge was stored reverse in (vertex_low, vertex_high).
+    int  face_edge;
+    bool operator==(const EdgeToFace &other) const { return vertex_low == other.vertex_low && vertex_high == other.vertex_high; }
+    bool operator<(const EdgeToFace &other) const { return vertex_low < other.vertex_low || (vertex_low == other.vertex_low && vertex_high < other.vertex_high); }
+};
 
-    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;
-}
-
-void VertexFaceIndex::create(const indexed_triangle_set &its)
-{
-    m_vertex_to_face_start.assign(its.vertices.size() + 1, 0);
-    // 1) Calculate vertex incidence by scatter.
-    for (auto &face : its.indices) {
-        ++ m_vertex_to_face_start[face(0) + 1];
-        ++ m_vertex_to_face_start[face(1) + 1];
-        ++ m_vertex_to_face_start[face(2) + 1];
-    }
-    // 2) Prefix sum to calculate offsets to m_vertex_faces_all.
-    for (size_t i = 2; i < m_vertex_to_face_start.size(); ++ i)
-        m_vertex_to_face_start[i] += m_vertex_to_face_start[i - 1];
-    // 3) Scatter indices of faces incident to a vertex into m_vertex_faces_all.
-    m_vertex_faces_all.assign(m_vertex_to_face_start.back(), 0);
-    for (size_t face_idx = 0; face_idx < its.indices.size(); ++ face_idx) {
-        auto &face = its.indices[face_idx];
-        for (int i = 0; i < 3; ++ i)
-            m_vertex_faces_all[m_vertex_to_face_start[face(i)] ++] = face_idx;
-    }
-    // 4) The previous loop modified m_vertex_to_face_start. Revert the change.
-    for (auto i = int(m_vertex_to_face_start.size()) - 1; i > 0; -- i)
-        m_vertex_to_face_start[i] = m_vertex_to_face_start[i - 1];
-    m_vertex_to_face_start.front() = 0;
-}
-
-// 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<Vec3i> create_face_neighbors_index_impl(const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
+static std::vector<EdgeToFace> create_edge_map(
+    const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
 {
-    std::vector<Vec3i> out(its.indices.size(), Vec3i(-1, -1, -1));
-
-    // Create a mapping from triangle edge into face.
-    struct EdgeToFace {
-        // Index of the 1st vertex of the triangle edge. vertex_low <= vertex_high.
-        int  vertex_low;
-        // Index of the 2nd vertex of the triangle edge.
-        int  vertex_high;
-        // Index of a triangular face.
-        int  face;
-        // Index of edge in the face, starting with 1. Negative indices if the edge was stored reverse in (vertex_low, vertex_high).
-        int  face_edge;
-        bool operator==(const EdgeToFace &other) const { return vertex_low == other.vertex_low && vertex_high == other.vertex_high; }
-        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(its.indices.size() * 3, EdgeToFace());
     for (uint32_t facet_idx = 0; facet_idx < its.indices.size(); ++ facet_idx)
@@ -750,6 +722,18 @@ static inline std::vector<Vec3i> create_face_neighbors_index_impl(const indexed_
     throw_on_cancel();
     std::sort(edges_map.begin(), edges_map.end());
 
+    return edges_map;
+}
+
+// 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<Vec3i> create_face_neighbors_index_impl(const indexed_triangle_set &its, ThrowOnCancelCallback throw_on_cancel)
+{
+    std::vector<Vec3i> out(its.indices.size(), Vec3i(-1, -1, -1));
+
+    std::vector<EdgeToFace> edges_map = create_edge_map(its, throw_on_cancel);
+
     // Assign a unique common edge id to touching triangle edges.
     int num_edges = 0;
     for (size_t i = 0; i < edges_map.size(); ++ i) {
@@ -981,20 +965,22 @@ Polygon its_convex_hull_2d_above(const indexed_triangle_set &its, const Transfor
 }
 
 // Generate the vertex list for a cube solid of arbitrary size in X/Y/Z.
+indexed_triangle_set its_make_cube(double xd, double yd, double zd)
+{
+    auto x = float(xd), y = float(yd), z = float(zd);
+    indexed_triangle_set mesh;
+    mesh.vertices = {{x, y, 0}, {x, 0, 0}, {0, 0, 0}, {0, y, 0},
+                     {x, y, z}, {0, y, z}, {0, 0, z}, {x, 0, z}};
+    mesh.indices  = {{0, 1, 2}, {0, 2, 3}, {4, 5, 6}, {4, 6, 7},
+                    {0, 4, 7}, {0, 7, 1}, {1, 7, 6}, {1, 6, 2},
+                    {2, 6, 5}, {2, 5, 3}, {4, 0, 3}, {4, 3, 5}};
+
+    return mesh;
+}
+
 TriangleMesh make_cube(double x, double y, double z) 
 {
-    TriangleMesh mesh(
-        {
-            {x, y, 0}, {x, 0, 0}, {0, 0, 0},
-            {0, y, 0}, {x, y, z}, {0, y, z},
-            {0, 0, z}, {x, 0, z}
-        },
-        {
-            {0, 1, 2}, {0, 2, 3}, {4, 5, 6},
-            {4, 6, 7}, {0, 4, 7}, {0, 7, 1},
-            {1, 7, 6}, {1, 6, 2}, {2, 6, 5},
-            {2, 5, 3}, {4, 0, 3}, {4, 3, 5}
-        });
+    TriangleMesh mesh(its_make_cube(x, y, z));
     mesh.repair();
     return mesh;
 }
@@ -1002,31 +988,32 @@ TriangleMesh make_cube(double x, double y, double z)
 // Generate the mesh for a cylinder and return it, using 
 // the generated angle to calculate the top mesh triangles.
 // Default is 360 sides, angle fa is in radians.
-TriangleMesh make_cylinder(double r, double h, double fa)
+indexed_triangle_set its_make_cylinder(double r, double h, double fa)
 {
+    indexed_triangle_set mesh;
     size_t n_steps    = (size_t)ceil(2. * PI / fa);
     double angle_step = 2. * PI / n_steps;
 
-    Pointf3s            vertices;
-    std::vector<Vec3i>  facets;
+    auto &vertices = mesh.vertices;
+    auto &facets   = mesh.indices;
     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));
-    vertices.emplace_back(Vec3d(0.0, 0.0, h));
+    vertices.emplace_back(Vec3f(0.f, 0.f, 0.f));
+    vertices.emplace_back(Vec3f(0.f, 0.f, float(h)));
 
     // for each line along the polygon approximating the top/bottom of the
     // 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));
+    Vec2f p = Eigen::Rotation2Df(0.f) * Eigen::Vector2f(0, r);
+    vertices.emplace_back(Vec3f(p(0), p(1), 0.f));
+    vertices.emplace_back(Vec3f(p(0), p(1), float(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));
+        p = Eigen::Rotation2Df(angle_step * i) * Eigen::Vector2f(0, float(r));
+        vertices.emplace_back(Vec3f(p(0), p(1), 0.f));
+        vertices.emplace_back(Vec3f(p(0), p(1), float(h)));
         int id = (int)vertices.size() - 1;
         facets.emplace_back( 0, id - 1, id - 3); // top
         facets.emplace_back(id,      1, id - 2); // bottom
@@ -1039,9 +1026,15 @@ TriangleMesh make_cylinder(double r, double h, double fa)
     facets.emplace_back( 3, 1,     id);
     facets.emplace_back(id, 2,      3);
     facets.emplace_back(id, id - 1, 2);
-    
-    TriangleMesh mesh(std::move(vertices), std::move(facets));
+
+    return mesh;
+}
+
+TriangleMesh make_cylinder(double r, double h, double fa)
+{
+    TriangleMesh mesh{its_make_cylinder(r, h, fa)};
     mesh.repair();
+
     return mesh;
 }
 
@@ -1076,14 +1069,15 @@ TriangleMesh make_cone(double r, double h, double fa)
 // to determine the granularity. 
 // Default angle is 1 degree.
 //FIXME better to discretize an Icosahedron recursively http://www.songho.ca/opengl/gl_sphere.html
-TriangleMesh make_sphere(double radius, double fa)
+indexed_triangle_set its_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);
 
-    Pointf3s vertices;
+    indexed_triangle_set mesh;
+    auto& vertices = mesh.vertices;
     vertices.reserve((stackCount - 1) * sectorCount + 2);
     for (int i = 0; i <= stackCount; ++ i) {
         // from pi/2 to -pi/2
@@ -1091,16 +1085,16 @@ TriangleMesh make_sphere(double radius, double fa)
         double xy = radius * cos(stackAngle);
         double z  = radius * sin(stackAngle);
         if (i == 0 || i == stackCount)
-            vertices.emplace_back(Vec3d(xy, 0., z));
+            vertices.emplace_back(Vec3f(float(xy), 0.f, float(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));
+                vertices.emplace_back(Vec3d(xy * std::cos(sectorAngle), xy * std::sin(sectorAngle), z).cast<float>());
             }
     }
 
-    std::vector<Vec3i> facets;
+    auto& facets = mesh.indices;
     facets.reserve(2 * (stackCount - 1) * sectorCount);
     for (int i = 0; i < stackCount; ++ i) {
         // Beginning of current stack.
@@ -1125,9 +1119,114 @@ TriangleMesh make_sphere(double radius, double fa)
             k2 = k2_next;
         }
     }
-    TriangleMesh mesh(std::move(vertices), std::move(facets));
-    mesh.repair();
+
     return mesh;
 }
 
+TriangleMesh make_sphere(double radius, double fa)
+{
+    TriangleMesh mesh(its_make_sphere(radius, fa));
+    mesh.repair();
+
+    return mesh;
 }
+
+void its_merge(indexed_triangle_set &A, const indexed_triangle_set &B)
+{
+    auto N   = int(A.vertices.size());
+    auto N_f = A.indices.size();
+
+    A.vertices.insert(A.vertices.end(), B.vertices.begin(), B.vertices.end());
+    A.indices.insert(A.indices.end(), B.indices.begin(), B.indices.end());
+
+    for(size_t n = N_f; n < A.indices.size(); n++)
+        A.indices[n] += Vec3i{N, N, N};
+}
+
+void its_merge(indexed_triangle_set &A, const std::vector<Vec3f> &triangles)
+{
+    const size_t offs = A.vertices.size();
+    A.vertices.insert(A.vertices.end(), triangles.begin(), triangles.end());
+    A.indices.reserve(A.indices.size() + A.vertices.size() / 3);
+
+    for(int i = int(offs); i < int(A.vertices.size()); i += 3)
+        A.indices.emplace_back(i, i + 1, i + 2);
+}
+
+void its_merge(indexed_triangle_set &A, const Pointf3s &triangles)
+{
+    auto trianglesf = reserve_vector<Vec3f> (triangles.size());
+    for (auto &t : triangles)
+        trianglesf.emplace_back(t.cast<float>());
+
+    its_merge(A, trianglesf);
+}
+
+float its_volume(const indexed_triangle_set &its)
+{
+    if (its.empty()) return 0.;
+
+    // Choose a point, any point as the reference.
+    auto p0 = its.vertices.front();
+    float volume = 0.f;
+    for (size_t i = 0; i < its.indices.size(); ++ i) {
+        // Do dot product to get distance from point to plane.
+        its_triangle triangle = its_triangle_vertices(its, i);
+        Vec3f U = triangle[1] - triangle[0];
+        Vec3f V = triangle[2] - triangle[0];
+        Vec3f C = U.cross(V);
+        Vec3f normal = C.normalized();
+        float area = 0.5 * C.norm();
+        float height = normal.dot(triangle[0] - p0);
+        volume += (area * height) / 3.0f;
+    }
+
+    return volume;
+}
+
+std::vector<indexed_triangle_set> its_split(const indexed_triangle_set &its)
+{
+    return its_split<>(its);
+}
+
+bool its_is_splittable(const indexed_triangle_set &its)
+{
+    return its_is_splittable<>(its);
+}
+
+void VertexFaceIndex::create(const indexed_triangle_set &its)
+{
+    m_vertex_to_face_start.assign(its.vertices.size() + 1, 0);
+    // 1) Calculate vertex incidence by scatter.
+    for (auto &face : its.indices) {
+        ++ m_vertex_to_face_start[face(0) + 1];
+        ++ m_vertex_to_face_start[face(1) + 1];
+        ++ m_vertex_to_face_start[face(2) + 1];
+    }
+    // 2) Prefix sum to calculate offsets to m_vertex_faces_all.
+    for (size_t i = 2; i < m_vertex_to_face_start.size(); ++ i)
+        m_vertex_to_face_start[i] += m_vertex_to_face_start[i - 1];
+    // 3) Scatter indices of faces incident to a vertex into m_vertex_faces_all.
+    m_vertex_faces_all.assign(m_vertex_to_face_start.back(), 0);
+    for (size_t face_idx = 0; face_idx < its.indices.size(); ++ face_idx) {
+        auto &face = its.indices[face_idx];
+        for (int i = 0; i < 3; ++ i)
+            m_vertex_faces_all[m_vertex_to_face_start[face(i)] ++] = face_idx;
+    }
+    // 4) The previous loop modified m_vertex_to_face_start. Revert the change.
+    for (auto i = int(m_vertex_to_face_start.size()) - 1; i > 0; -- i)
+        m_vertex_to_face_start[i] = m_vertex_to_face_start[i - 1];
+    m_vertex_to_face_start.front() = 0;
+}
+
+std::vector<Vec3i> its_create_neighbors_index(const indexed_triangle_set &its)
+{
+    return create_neighbors_index(ex_seq, its);
+}
+
+std::vector<Vec3i> its_create_neighbors_index_par(const indexed_triangle_set &its)
+{
+    return create_neighbors_index(ex_tbb, its);
+}
+
+} // namespace Slic3r
diff --git a/src/libslic3r/TriangleMesh.hpp b/src/libslic3r/TriangleMesh.hpp
index 6d7fd48e4f..8ac0e98717 100644
--- a/src/libslic3r/TriangleMesh.hpp
+++ b/src/libslic3r/TriangleMesh.hpp
@@ -88,11 +88,6 @@ private:
     std::deque<uint32_t> find_unvisited_neighbors(std::vector<unsigned char> &facet_visited) const;
 };
 
-// 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_vertex_faces_index(const indexed_triangle_set &its);
-
 // Index of face indices incident with a vertex index.
 struct VertexFaceIndex
 {
@@ -111,6 +106,8 @@ public:
     // Vertex incidence.
     size_t   count(size_t vertex_id) const throw() { return m_vertex_to_face_start[vertex_id + 1] - m_vertex_to_face_start[vertex_id]; }
 
+    const Range<iterator> operator[](size_t vertex_id) const { return {begin(vertex_id), end(vertex_id)}; }
+
 private:
     std::vector<size_t>     m_vertex_to_face_start;
     std::vector<size_t>     m_vertex_faces_all;
@@ -122,6 +119,11 @@ private:
 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);
 
+// Create index that gives neighbor faces for each face. Ignores face orientations.
+// TODO: naming...
+std::vector<Vec3i> its_create_neighbors_index(const indexed_triangle_set &its);
+std::vector<Vec3i> its_create_neighbors_index_par(const indexed_triangle_set &its);
+
 // After applying a transformation with negative determinant, flip the faces to keep the transformed mesh volume positive.
 void its_flip_triangles(indexed_triangle_set &its);
 
@@ -136,6 +138,10 @@ int its_remove_degenerate_faces(indexed_triangle_set &its, bool shrink_to_fit =
 // 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);
 
+std::vector<indexed_triangle_set> its_split(const indexed_triangle_set &its);
+
+bool its_is_splittable(const indexed_triangle_set &its);
+
 // Shrink the vectors of its.vertices and its.faces to a minimum size by reallocating the two vectors.
 void its_shrink_to_fit(indexed_triangle_set &its);
 
@@ -147,11 +153,56 @@ void its_collect_mesh_projection_points_above(const indexed_triangle_set &its, c
 Polygon its_convex_hull_2d_above(const indexed_triangle_set &its, const Matrix3f &m, const float z);
 Polygon its_convex_hull_2d_above(const indexed_triangle_set &its, const Transform3f &t, const float z);
 
+using its_triangle = std::array<stl_vertex, 3>;
+
+inline its_triangle its_triangle_vertices(const indexed_triangle_set &its,
+                                          size_t                      face_id)
+{
+    return {its.vertices[its.indices[face_id](0)],
+            its.vertices[its.indices[face_id](1)],
+            its.vertices[its.indices[face_id](2)]};
+}
+
+inline stl_normal its_unnormalized_normal(const indexed_triangle_set &its,
+                                          size_t                      face_id)
+{
+    its_triangle tri = its_triangle_vertices(its, face_id);
+    return (tri[1] - tri[0]).cross(tri[2] - tri[0]);
+}
+
+float its_volume(const indexed_triangle_set &its);
+
+void its_merge(indexed_triangle_set &A, const indexed_triangle_set &B);
+void its_merge(indexed_triangle_set &A, const std::vector<Vec3f> &triangles);
+void its_merge(indexed_triangle_set &A, const Pointf3s &triangles);
+
+indexed_triangle_set its_make_cube(double x, double y, double z);
 TriangleMesh make_cube(double x, double y, double z);
+
+// Generate a TriangleMesh of a cylinder
+indexed_triangle_set its_make_cylinder(double r, double h, double fa=(2*PI/360));
 TriangleMesh make_cylinder(double r, double h, double fa=(2*PI/360));
+
+indexed_triangle_set its_make_sphere(double rho, double fa=(2*PI/360));
 TriangleMesh make_cone(double r, double h, double fa=(2*PI/360));
 TriangleMesh make_sphere(double rho, double fa=(2*PI/360));
 
+inline BoundingBoxf3 bounding_box(const TriangleMesh &m) { return m.bounding_box(); }
+inline BoundingBoxf3 bounding_box(const indexed_triangle_set& its)
+{
+    if (its.vertices.empty())
+        return {};
+
+    Vec3f bmin = its.vertices.front(), bmax = its.vertices.front();
+
+    for (const Vec3f &p : its.vertices) {
+        bmin = p.cwiseMin(bmin);
+        bmax = p.cwiseMax(bmax);
+    }
+
+    return {bmin.cast<double>(), bmax.cast<double>()};
+}
+
 }
 
 // Serialization through the Cereal library
diff --git a/src/libslic3r/TriangulateWall.cpp b/src/libslic3r/TriangulateWall.cpp
index ec2945b10b..eb710d9a57 100644
--- a/src/libslic3r/TriangulateWall.cpp
+++ b/src/libslic3r/TriangulateWall.cpp
@@ -3,131 +3,157 @@
 
 namespace Slic3r {
 
-class Ring {
-    size_t idx = 0, nextidx = 1, startidx = 0, begin = 0, end = 0;
+//class Ring {
+//    size_t idx = 0, nextidx = 1, startidx = 0, begin = 0, end = 0;
     
-public:
-    explicit Ring(size_t from, size_t to) : begin(from), end(to) { init(begin); }
+//public:
+//    explicit Ring(size_t from, size_t to) : begin(from), end(to) { init(begin); }
 
-    size_t size() const { return end - begin; }
-    std::pair<size_t, size_t> pos() const { return {idx, nextidx}; }
-    bool is_lower() const { return idx < size(); }
+//    size_t size() const { return end - begin; }
+//    std::pair<size_t, size_t> pos() const { return {idx, nextidx}; }
+//    bool is_lower() const { return idx < size(); }
     
-    void inc()
-    {
-        if (nextidx != startidx) nextidx++;
-        if (nextidx == end) nextidx = begin;
-        idx ++;
-        if (idx == end) idx = begin;
-    }
+//    void inc()
+//    {
+//        if (nextidx != startidx) nextidx++;
+//        if (nextidx == end) nextidx = begin;
+//        idx ++;
+//        if (idx == end) idx = begin;
+//    }
     
-    void init(size_t pos)
-    {
-        startidx = begin + (pos - begin) % size();
-        idx = startidx;
-        nextidx = begin + (idx + 1 - begin) % size();
-    }
+//    void init(size_t pos)
+//    {
+//        startidx = begin + (pos - begin) % size();
+//        idx = startidx;
+//        nextidx = begin + (idx + 1 - begin) % size();
+//    }
     
-    bool is_finished() const { return nextidx == idx; }
-};
+//    bool is_finished() const { return nextidx == idx; }
+//};
 
-static double sq_dst(const Vec3d &v1, const Vec3d& v2)
-{
-    Vec3d v = v1 - v2;
-    return v.x() * v.x() + v.y() * v.y() /*+ v.z() * v.z()*/;
-}
+//template<class Sc>
+//static Sc sq_dst(const Vec<3, Sc> &v1, const Vec<3, Sc>& v2)
+//{
+//    Vec<3, Sc> v = v1 - v2;
+//    return v.x() * v.x() + v.y() * v.y() /*+ v.z() * v.z()*/;
+//}
 
-static double score(const Ring& onring, const Ring &offring,
-                    const std::vector<Vec3d> &pts)
-{
-    double a = sq_dst(pts[onring.pos().first], pts[offring.pos().first]);
-    double b = sq_dst(pts[onring.pos().second], pts[offring.pos().first]);
-    return (std::abs(a) + std::abs(b)) / 2.;
-}
+//template<class Sc>
+//static Sc trscore(const Ring &                   onring,
+//                  const Ring &                   offring,
+//                  const std::vector<Vec<3, Sc>> &pts)
+//{
+//    Sc a = sq_dst(pts[onring.pos().first], pts[offring.pos().first]);
+//    Sc b = sq_dst(pts[onring.pos().second], pts[offring.pos().first]);
+//    return (std::abs(a) + std::abs(b)) / 2.;
+//}
 
-class Triangulator {
-    const std::vector<Vec3d> *pts;
-    Ring *onring, *offring;
+//template<class Sc>
+//class Triangulator {
+//    const std::vector<Vec<3, Sc>> *pts;
+//    Ring *onring, *offring;
     
-    double calc_score() const
-    {
-        return Slic3r::score(*onring, *offring, *pts);
-    }
+//    double calc_score() const
+//    {
+//        return trscore(*onring, *offring, *pts);
+//    }
     
-    void synchronize_rings()
-    {
-        Ring lring = *offring;
-        auto minsc = Slic3r::score(*onring, lring, *pts);
-        size_t imin = lring.pos().first;
+//    void synchronize_rings()
+//    {
+//        Ring lring = *offring;
+//        auto minsc = trscore(*onring, lring, *pts);
+//        size_t imin = lring.pos().first;
         
-        lring.inc();
+//        lring.inc();
         
-        while(!lring.is_finished()) {
-            double score = Slic3r::score(*onring, lring, *pts);
-            if (score < minsc) { minsc = score; imin = lring.pos().first; }
-            lring.inc();
-        }
+//        while(!lring.is_finished()) {
+//            double score = trscore(*onring, lring, *pts);
+//            if (score < minsc) { minsc = score; imin = lring.pos().first; }
+//            lring.inc();
+//        }
         
-        offring->init(imin);
-    }
+//        offring->init(imin);
+//    }
     
-    void emplace_indices(std::vector<Vec3i> &indices)
-    {
-        Vec3i tr{int(onring->pos().first), int(onring->pos().second),
-                 int(offring->pos().first)};
-        if (onring->is_lower()) std::swap(tr(0), tr(1));
-        indices.emplace_back(tr);
-    }
+//    void emplace_indices(std::vector<Vec3i> &indices)
+//    {
+//        Vec3i tr{int(onring->pos().first), int(onring->pos().second),
+//                 int(offring->pos().first)};
+//        if (onring->is_lower()) std::swap(tr(0), tr(1));
+//        indices.emplace_back(tr);
+//    }
     
-public:
-    void run(std::vector<Vec3i> &indices)
-    {   
-        synchronize_rings();
+//public:
+//    void run(std::vector<Vec3i> &indices)
+//    {
+//        synchronize_rings();
         
-        double score = 0, prev_score = 0;        
-        while (!onring->is_finished() || !offring->is_finished()) {
-            prev_score = score;
-            if (onring->is_finished() || (score = calc_score()) > prev_score) {
-                std::swap(onring, offring);
-            } else {
-                emplace_indices(indices);
-                onring->inc();
-            }
-        }
-    }
+//        double score = 0, prev_score = 0;
+//        while (!onring->is_finished() || !offring->is_finished()) {
+//            prev_score = score;
+//            if (onring->is_finished() || (score = calc_score()) > prev_score) {
+//                std::swap(onring, offring);
+//            } else {
+//                emplace_indices(indices);
+//                onring->inc();
+//            }
+//        }
+//    }
 
-    explicit Triangulator(const std::vector<Vec3d> *points,
-                          Ring &                    lower,
-                          Ring &                    upper)
-        : pts{points}, onring{&upper}, offring{&lower}
-    {}
-};
+//    explicit Triangulator(const std::vector<Vec<3, Sc>> *points,
+//                          Ring &                    lower,
+//                          Ring &                    upper)
+//        : pts{points}, onring{&upper}, offring{&lower}
+//    {}
+//};
 
-Wall triangulate_wall(
-    const Polygon &          lower,
-    const Polygon &          upper,
-    double                   lower_z_mm,
-    double                   upper_z_mm)
-{
-    if (upper.points.size() < 3 || lower.points.size() < 3) return {};
-    
-    Wall wall;
-    auto &pts = wall.first;
-    auto &ind = wall.second;
+//template<class Sc, class I>
+//void triangulate_wall(std::vector<Vec<3, Sc>> &pts,
+//                      std::vector<Vec<3, I>> & ind,
+//                      const Polygon &          lower,
+//                      const Polygon &          upper,
+//                      double                   lower_z_mm,
+//                      double                   upper_z_mm)
+//{
+//    if (upper.points.size() < 3 || lower.points.size() < 3) return;
 
-    pts.reserve(lower.points.size() + upper.points.size());
-    for (auto &p : lower.points)
-        wall.first.emplace_back(unscaled(p.x()), unscaled(p.y()), lower_z_mm);
-    for (auto &p : upper.points)
-        wall.first.emplace_back(unscaled(p.x()), unscaled(p.y()), upper_z_mm);
+//    pts.reserve(lower.points.size() + upper.points.size());
+//    for (auto &p : lower.points)
+//        pts.emplace_back(unscaled(p.x()), unscaled(p.y()), lower_z_mm);
+//    for (auto &p : upper.points)
+//        pts.emplace_back(unscaled(p.x()), unscaled(p.y()), upper_z_mm);
+
+//    ind.reserve(2 * (lower.size() + upper.size()));
+
+//    Ring lring{0, lower.points.size()}, uring{lower.points.size(), pts.size()};
+//    Triangulator t{&pts, lring, uring};
+//    t.run(ind);
+//}
+
+//Wall triangulate_wall(const Polygon &lower,
+//                      const Polygon &upper,
+//                      double         lower_z_mm,
+//                      double         upper_z_mm)
+//{
+//    if (upper.points.size() < 3 || lower.points.size() < 3) return {};
     
-    ind.reserve(2 * (lower.size() + upper.size()));
+//    Wall wall;
+//    auto &pts = wall.first;
+//    auto &ind = wall.second;
+
+//    pts.reserve(lower.points.size() + upper.points.size());
+//    for (auto &p : lower.points)
+//        wall.first.emplace_back(unscaled(p.x()), unscaled(p.y()), lower_z_mm);
+//    for (auto &p : upper.points)
+//        wall.first.emplace_back(unscaled(p.x()), unscaled(p.y()), upper_z_mm);
     
-    Ring lring{0, lower.points.size()}, uring{lower.points.size(), pts.size()};
-    Triangulator t{&pts, lring, uring};
-    t.run(ind);
+//    ind.reserve(2 * (lower.size() + upper.size()));
     
-    return wall;
-}
+//    Ring lring{0, lower.points.size()}, uring{lower.points.size(), pts.size()};
+//    Triangulator t{&pts, lring, uring};
+//    t.run(ind);
+    
+//    return wall;
+//}
 
 } // namespace Slic3r
diff --git a/src/libslic3r/TriangulateWall.hpp b/src/libslic3r/TriangulateWall.hpp
index 68bf4b0ac6..c64a75e63c 100644
--- a/src/libslic3r/TriangulateWall.hpp
+++ b/src/libslic3r/TriangulateWall.hpp
@@ -5,13 +5,148 @@
 
 namespace Slic3r {
 
-using Wall = std::pair<std::vector<Vec3d>, std::vector<Vec3i>>;
+namespace trianglulate_wall_detail {
 
-Wall triangulate_wall(
-    const Polygon &       lower,
-    const Polygon &       upper,
-    double                lower_z_mm,
-    double                upper_z_mm);
+class Ring {
+    size_t idx = 0, nextidx = 1, startidx = 0, begin = 0, end = 0;
+
+public:
+    explicit Ring(size_t from, size_t to) : begin(from), end(to) { init(begin); }
+
+    size_t size() const { return end - begin; }
+    std::pair<size_t, size_t> pos() const { return {idx, nextidx}; }
+    bool is_lower() const { return idx < size(); }
+
+    void inc()
+    {
+        if (nextidx != startidx) nextidx++;
+        if (nextidx == end) nextidx = begin;
+        idx ++;
+        if (idx == end) idx = begin;
+    }
+
+    void init(size_t pos)
+    {
+        startidx = begin + (pos - begin) % size();
+        idx = startidx;
+        nextidx = begin + (idx + 1 - begin) % size();
+    }
+
+    bool is_finished() const { return nextidx == idx; }
+};
+
+template<class Sc>
+static Sc sq_dst(const Vec<3, Sc> &v1, const Vec<3, Sc>& v2)
+{
+    Vec<3, Sc> v = v1 - v2;
+    return v.x() * v.x() + v.y() * v.y() /*+ v.z() * v.z()*/;
 }
 
+template<class Sc>
+static Sc trscore(const Ring &                   onring,
+                  const Ring &                   offring,
+                  const std::vector<Vec<3, Sc>> &pts)
+{
+    Sc a = sq_dst(pts[onring.pos().first], pts[offring.pos().first]);
+    Sc b = sq_dst(pts[onring.pos().second], pts[offring.pos().first]);
+    return (std::abs(a) + std::abs(b)) / 2.;
+}
+
+template<class Sc>
+class Triangulator {
+    const std::vector<Vec<3, Sc>> *pts;
+    Ring *onring, *offring;
+
+    double calc_score() const
+    {
+        return trscore(*onring, *offring, *pts);
+    }
+
+    void synchronize_rings()
+    {
+        Ring lring = *offring;
+        auto minsc = trscore(*onring, lring, *pts);
+        size_t imin = lring.pos().first;
+
+        lring.inc();
+
+        while(!lring.is_finished()) {
+            double score = trscore(*onring, lring, *pts);
+            if (score < minsc) { minsc = score; imin = lring.pos().first; }
+            lring.inc();
+        }
+
+        offring->init(imin);
+    }
+
+    void emplace_indices(std::vector<Vec3i> &indices)
+    {
+        Vec3i tr{int(onring->pos().first), int(onring->pos().second),
+                 int(offring->pos().first)};
+        if (onring->is_lower()) std::swap(tr(0), tr(1));
+        indices.emplace_back(tr);
+    }
+
+public:
+    void run(std::vector<Vec3i> &indices)
+    {
+        synchronize_rings();
+
+        double score = 0, prev_score = 0;
+        while (!onring->is_finished() || !offring->is_finished()) {
+            prev_score = score;
+            if (onring->is_finished() || (score = calc_score()) > prev_score) {
+                std::swap(onring, offring);
+            } else {
+                emplace_indices(indices);
+                onring->inc();
+            }
+        }
+    }
+
+    explicit Triangulator(const std::vector<Vec<3, Sc>> *points,
+                          Ring &                    lower,
+                          Ring &                    upper)
+        : pts{points}, onring{&upper}, offring{&lower}
+    {}
+};
+
+} // namespace trianglulate_wall_detail
+
+template<class Sc, class I>
+void triangulate_wall(std::vector<Vec<3, Sc>> &pts,
+                      std::vector<Vec<3, I>> & ind,
+                      const Polygon &          lower,
+                      const Polygon &          upper,
+                      double                   lower_z_mm,
+                      double                   upper_z_mm)
+{
+    using namespace trianglulate_wall_detail;
+
+    if (upper.points.size() < 3 || lower.points.size() < 3) return;
+
+    pts.reserve(lower.points.size() + upper.points.size());
+    for (auto &p : lower.points)
+        pts.emplace_back(unscaled(p.x()), unscaled(p.y()), lower_z_mm);
+    for (auto &p : upper.points)
+        pts.emplace_back(unscaled(p.x()), unscaled(p.y()), upper_z_mm);
+
+    ind.reserve(2 * (lower.size() + upper.size()));
+
+    Ring lring{0, lower.points.size()}, uring{lower.points.size(), pts.size()};
+    Triangulator t{&pts, lring, uring};
+    t.run(ind);
+}
+
+//using Wall = std::pair<std::vector<Vec3d>, std::vector<Vec3i>>;
+
+//Wall triangulate_wall(
+//    const Polygon &       lower,
+//    const Polygon &       upper,
+//    double                lower_z_mm,
+//    double                upper_z_mm);
+//}
+
+} // namespace Slic3r
+
 #endif // TRIANGULATEWALL_HPP
diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h
index d2c1d0fc18..8d8d15df9d 100644
--- a/src/libslic3r/libslic3r.h
+++ b/src/libslic3r/libslic3r.h
@@ -306,6 +306,29 @@ IntegerOnly<I, std::vector<T, Args...>> reserve_vector(I capacity)
 template<class T>
 using remove_cvref_t = std::remove_cv_t<std::remove_reference_t<T>>;
 
+// A very simple range concept implementation with iterator-like objects.
+// This should be replaced by std::ranges::subrange (C++20)
+template<class It> class Range
+{
+    It from, to;
+public:
+
+    // The class is ready for range based for loops.
+    It begin() const { return from; }
+    It end() const { return to; }
+
+    // The iterator type can be obtained this way.
+    using iterator = It;
+    using value_type = typename std::iterator_traits<It>::value_type;
+
+    Range() = default;
+    Range(It b, It e) : from(std::move(b)), to(std::move(e)) {}
+
+    // Some useful container-like methods...
+    inline size_t size() const { return end() - begin(); }
+    inline bool   empty() const { return size() == 0; }
+};
+
 } // namespace Slic3r
 
 #endif
diff --git a/src/slic3r/GUI/Gizmos/GLGizmosCommon.cpp b/src/slic3r/GUI/Gizmos/GLGizmosCommon.cpp
index acb85d5399..425e11f73e 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmosCommon.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmosCommon.cpp
@@ -208,7 +208,7 @@ void HollowedMesh::on_update()
                     m_drainholes = print_object->model_object()->sla_drain_holes;
                     m_old_hollowing_timestamp = timestamp;
 
-                    const TriangleMesh &interior = print_object->hollowed_interior_mesh();
+                    const indexed_triangle_set &interior = print_object->hollowed_interior_mesh();
                     if (!interior.empty()) {
                         m_hollowed_interior_transformed = std::make_unique<TriangleMesh>(interior);
                         m_hollowed_interior_transformed->repaired = false;
diff --git a/src/slic3r/GUI/Jobs/SLAImportJob.cpp b/src/slic3r/GUI/Jobs/SLAImportJob.cpp
index 3da937f0ac..9998f42c76 100644
--- a/src/slic3r/GUI/Jobs/SLAImportJob.cpp
+++ b/src/slic3r/GUI/Jobs/SLAImportJob.cpp
@@ -113,14 +113,14 @@ public:
     Plater *plater;
     
     Sel sel = Sel::modelAndProfile;
-    
-    TriangleMesh       mesh;
-    DynamicPrintConfig profile;
-    wxString           path;
-    Vec2i              win = {2, 2};
-    std::string        err;
-    
-    priv(Plater *plt): plater{plt} {}
+
+    indexed_triangle_set mesh;
+    DynamicPrintConfig   profile;
+    wxString             path;
+    Vec2i                win = {2, 2};
+    std::string          err;
+
+    priv(Plater *plt) : plater{plt} {}
 };
 
 SLAImportJob::SLAImportJob(std::shared_ptr<ProgressIndicator> pri, Plater *plater)
@@ -222,7 +222,8 @@ void SLAImportJob::finalize()
     
     if (!p->mesh.empty()) {
         bool is_centered = false;
-        p->plater->sidebar().obj_list()->load_mesh_object(p->mesh, name, is_centered);
+        p->plater->sidebar().obj_list()->load_mesh_object(TriangleMesh{p->mesh},
+                                                          name, is_centered);
     }
     
     reset();
diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt
index f93cf10afd..5a4203ea0c 100644
--- a/tests/libslic3r/CMakeLists.txt
+++ b/tests/libslic3r/CMakeLists.txt
@@ -21,6 +21,7 @@ add_executable(${_TEST_NAME}_tests
     test_optimizers.cpp
     test_png_io.cpp
     test_timeutils.cpp
+    test_indexed_triangle_set.cpp
 	)
 
 if (TARGET OpenVDB::openvdb)
diff --git a/tests/libslic3r/test_indexed_triangle_set.cpp b/tests/libslic3r/test_indexed_triangle_set.cpp
new file mode 100644
index 0000000000..4c11283451
--- /dev/null
+++ b/tests/libslic3r/test_indexed_triangle_set.cpp
@@ -0,0 +1,102 @@
+#include <iostream>
+#include <fstream>
+#include <catch2/catch.hpp>
+
+#include "libslic3r/TriangleMesh.hpp"
+
+TEST_CASE("Split empty mesh", "[its_split][its]") {
+    using namespace Slic3r;
+
+    indexed_triangle_set its;
+
+    std::vector<indexed_triangle_set> res = its_split(its);
+
+    REQUIRE(res.empty());
+}
+
+TEST_CASE("Split simple mesh consisting of one part", "[its_split][its]") {
+    using namespace Slic3r;
+
+    auto cube = its_make_cube(10., 10., 10.);
+
+    std::vector<indexed_triangle_set> res = its_split(cube);
+
+    REQUIRE(res.size() == 1);
+    REQUIRE(res.front().indices.size() == cube.indices.size());
+    REQUIRE(res.front().vertices.size() == cube.vertices.size());
+}
+
+void debug_write_obj(const std::vector<indexed_triangle_set> &res, const std::string &name)
+{
+#ifndef NDEBUG
+    size_t part_idx = 0;
+    for (auto &part_its : res) {
+        its_write_obj(part_its, (name + std::to_string(part_idx++) + ".obj").c_str());
+    }
+#endif
+}
+
+TEST_CASE("Split two non-watertight mesh", "[its_split][its]") {
+    using namespace Slic3r;
+
+    auto cube1 = its_make_cube(10., 10., 10.);
+    cube1.indices.pop_back();
+    auto cube2 = cube1;
+
+    its_transform(cube1, identity3f().translate(Vec3f{-5.f, 0.f, 0.f}));
+    its_transform(cube2, identity3f().translate(Vec3f{5.f, 0.f, 0.f}));
+
+    its_merge(cube1, cube2);
+
+    std::vector<indexed_triangle_set> res = its_split(cube1);
+
+    REQUIRE(res.size() == 2);
+    REQUIRE(res[0].indices.size() == res[1].indices.size());
+    REQUIRE(res[0].indices.size() == cube2.indices.size());
+    REQUIRE(res[0].vertices.size() == res[1].vertices.size());
+    REQUIRE(res[0].vertices.size() == cube2.vertices.size());
+
+    debug_write_obj(res, "parts_non_watertight");
+}
+
+TEST_CASE("Split non-manifold mesh", "[its_split][its]") {
+    using namespace Slic3r;
+
+    auto cube = its_make_cube(10., 10., 10.), cube_low = cube;
+
+    its_transform(cube_low, identity3f().translate(Vec3f{10.f, 10.f, 10.f}));
+    its_merge(cube, cube_low);
+    its_merge_vertices(cube);
+
+    std::vector<indexed_triangle_set> res = its_split(cube);
+
+    REQUIRE(res.size() == 2);
+    REQUIRE(res[0].indices.size() == res[1].indices.size());
+    REQUIRE(res[0].indices.size() == cube_low.indices.size());
+    REQUIRE(res[0].vertices.size() == res[1].vertices.size());
+    REQUIRE(res[0].vertices.size() == cube_low.vertices.size());
+
+    debug_write_obj(res, "cubes_non_manifold");
+}
+
+TEST_CASE("Split two watertight meshes", "[its_split][its]") {
+    using namespace Slic3r;
+
+    auto sphere1 = its_make_sphere(10., 2 * PI / 200.), sphere2 = sphere1;
+
+    its_transform(sphere1, identity3f().translate(Vec3f{-5.f, 0.f, 0.f}));
+    its_transform(sphere2, identity3f().translate(Vec3f{5.f, 0.f, 0.f}));
+
+    its_merge(sphere1, sphere2);
+
+    std::vector<indexed_triangle_set> res = its_split(sphere1);
+
+    REQUIRE(res.size() == 2);
+    REQUIRE(res[0].indices.size() == res[1].indices.size());
+    REQUIRE(res[0].indices.size() == sphere2.indices.size());
+    REQUIRE(res[0].vertices.size() == res[1].vertices.size());
+    REQUIRE(res[0].vertices.size() == sphere2.vertices.size());
+
+    debug_write_obj(res, "parts_watertight");
+}
+
diff --git a/tests/libslic3r/test_marchingsquares.cpp b/tests/libslic3r/test_marchingsquares.cpp
index 14481f14a0..a0b27f6d87 100644
--- a/tests/libslic3r/test_marchingsquares.cpp
+++ b/tests/libslic3r/test_marchingsquares.cpp
@@ -17,7 +17,6 @@
 #include <libslic3r/TriangulateWall.hpp>
 #include <libslic3r/Tesselate.hpp>
 #include <libslic3r/SlicesToTriangleMesh.hpp>
-#include <libslic3r/SLA/Contour3D.hpp>
 
 using namespace Slic3r;
 
@@ -368,10 +367,9 @@ static void recreate_object_from_rasters(const std::string &objname, float lh) {
         layer = std::move(layer_);
     }
     
-    TriangleMesh out = slices_to_triangle_mesh(layers, bb.min.z(), double(lh), double(lh));
-    
-    out.require_shared_vertices();
-    out.WriteOBJFile("out_from_rasters.obj");
+    indexed_triangle_set out = slices_to_mesh(layers, bb.min.z(), double(lh), double(lh));
+
+    its_write_obj(out, "out_from_rasters.obj");
 }
 
 TEST_CASE("Recreate object from rasters", "[SL1Import]") {
diff --git a/tests/sla_print/sla_print_tests.cpp b/tests/sla_print/sla_print_tests.cpp
index 54e5ea7321..d26d6e2a23 100644
--- a/tests/sla_print/sla_print_tests.cpp
+++ b/tests/sla_print/sla_print_tests.cpp
@@ -223,23 +223,14 @@ TEST_CASE("RasterizedPolygonAreaShouldMatch", "[SLARasterOutput]") {
     REQUIRE(raster_pxsum(raster0) == 0);
 }
 
-TEST_CASE("Triangle mesh conversions should be correct", "[SLAConversions]")
-{
-    sla::Contour3D cntr;
-    
-    {
-        std::fstream infile{"extruder_idler_quads.obj", std::ios::in};
-        cntr.from_obj(infile);
-    }
-}
 
 TEST_CASE("halfcone test", "[halfcone]") {
     sla::DiffBridge br{Vec3d{1., 1., 1.}, Vec3d{10., 10., 10.}, 0.25, 0.5};
 
-    TriangleMesh m = sla::to_triangle_mesh(sla::get_mesh(br, 45));
+    indexed_triangle_set m = sla::get_mesh(br, 45);
 
-    m.require_shared_vertices();
-    m.WriteOBJFile("Halfcone.obj");
+    its_merge_vertices(m);
+    its_write_obj(m, "Halfcone.obj");
 }
 
 TEST_CASE("Test concurrency")
diff --git a/tests/sla_print/sla_test_utils.cpp b/tests/sla_print/sla_test_utils.cpp
index 8a78c63b3e..be9bf9741d 100644
--- a/tests/sla_print/sla_test_utils.cpp
+++ b/tests/sla_print/sla_test_utils.cpp
@@ -69,9 +69,10 @@ void export_failed_case(const std::vector<ExPolygons> &support_slices, const Sup
             svg.Close();
         }
     }
-    
-    TriangleMesh m;
-    byproducts.supporttree.retrieve_full_mesh(m);
+
+    indexed_triangle_set its;
+    byproducts.supporttree.retrieve_full_mesh(its);
+    TriangleMesh m{its};
     m.merge(byproducts.input_mesh);
     m.repair();
     m.require_shared_vertices();
@@ -93,7 +94,7 @@ void test_supports(const std::string          &obj_filename,
     if (hollowingcfg.enabled) {
         sla::InteriorPtr interior = sla::generate_interior(mesh, hollowingcfg);
         REQUIRE(interior);
-        mesh.merge(sla::get_mesh(*interior));
+        mesh.merge(TriangleMesh{sla::get_mesh(*interior)});
         mesh.require_shared_vertices();
     }
     
@@ -151,7 +152,7 @@ void test_supports(const std::string          &obj_filename,
     
     check_support_tree_integrity(treebuilder, supportcfg);
     
-    const TriangleMesh &output_mesh = treebuilder.retrieve_mesh();
+    TriangleMesh output_mesh{treebuilder.retrieve_mesh(sla::MeshType::Support)};
     
     check_validity(output_mesh, validityflags);
     
@@ -228,14 +229,16 @@ void test_pad(const std::string &obj_filename, const sla::PadConfig &padcfg, Pad
     REQUIRE_FALSE(mesh.empty());
     
     // Create pad skeleton only from the model
-    Slic3r::sla::pad_blueprint(mesh, out.model_contours);
+    Slic3r::sla::pad_blueprint(mesh.its, out.model_contours);
     
     test_concave_hull(out.model_contours);
     
     REQUIRE_FALSE(out.model_contours.empty());
     
     // Create the pad geometry for the model contours only
-    Slic3r::sla::create_pad({}, out.model_contours, out.mesh, padcfg);
+    indexed_triangle_set out_its;
+    Slic3r::sla::create_pad({}, out.model_contours, out_its, padcfg);
+    out.mesh = TriangleMesh{out_its};
     
     check_validity(out.mesh);