diff --git a/sandboxes/meshboolean/MeshBoolean.cpp b/sandboxes/meshboolean/MeshBoolean.cpp
index d339ef5c3..392d90707 100644
--- a/sandboxes/meshboolean/MeshBoolean.cpp
+++ b/sandboxes/meshboolean/MeshBoolean.cpp
@@ -11,66 +11,34 @@
 
 #include <boost/log/trivial.hpp>
 
-namespace Slic3r {
-
-} // namespace Slic3r
-
 int main(const int argc, const char * argv[])
 {
     using namespace Slic3r;
     
-    if (argc <= 1) return EXIT_FAILURE;
+    if (argc <= 1) {
+        std::cout << "Usage: meshboolean <input_file.3mf>" << std::endl;
+        return EXIT_FAILURE;
+    }
     
-    DynamicPrintConfig cfg;
-    auto model = Model::read_from_file(argv[1], &cfg);
     
-    if (model.objects.empty()) return EXIT_SUCCESS;
+    TriangleMesh input;
     
-    SLAPrint print;
-    print.apply(model, cfg);
-    PrintBase::TaskParams task;
-    task.to_object_step = slaposHollowing;
-    
-    print.set_task(task);
-    print.process();
+    input.ReadSTLFile(argv[1]);
+    input.repair();
     
     Benchmark bench;
     
-    for (SLAPrintObject *po : print.objects()) {
-        TriangleMesh holes;
-        sla::DrainHoles holepts = po->transformed_drainhole_points();
-        
-        for (auto &hole: holepts)
-            holes.merge(sla::to_triangle_mesh(hole.to_mesh()));
-        
-        TriangleMesh hollowed_mesh = po->transformed_mesh();
-        hollowed_mesh.merge(po->hollowed_interior_mesh());
-        
-        hollowed_mesh.require_shared_vertices();
-        holes.require_shared_vertices();
-        
-        TriangleMesh drilled_mesh_igl = hollowed_mesh;
-        bench.start();
-        MeshBoolean::minus(drilled_mesh_igl, holes);
-        bench.stop();
-        
-        std::cout << "Mesh boolean duration with IGL: " << bench.getElapsedSec() << std::endl;
-        
-        TriangleMesh drilled_mesh_cgal = hollowed_mesh;
-        bench.start();
-        MeshBoolean::cgal::self_union(drilled_mesh_cgal);
-        MeshBoolean::cgal::minus(drilled_mesh_cgal, holes);
-        bench.stop();
-        
-        std::cout << "Mesh boolean duration with CGAL: " << bench.getElapsedSec() << std::endl;
-        
-        std::string name("obj"), outf;
-        outf = name + "igl" + std::to_string(po->model_object()->id().id) + ".obj";
-        drilled_mesh_igl.WriteOBJFile(outf.c_str());
-        
-        outf = name + "cgal" + std::to_string(po->model_object()->id().id) + ".obj";
-        drilled_mesh_cgal.WriteOBJFile(outf.c_str());
-    }
+    bench.start();
+    bool fckd = MeshBoolean::cgal::does_self_intersect(input);
+    bench.stop();
+    
+    std::cout << "Self intersect test: " << fckd << " duration: " << bench.getElapsedSec() << std::endl;
+    
+    bench.start();
+    MeshBoolean::self_union(input);
+    bench.stop();
+    
+    std::cout << "Self union duration: " << bench.getElapsedSec() << std::endl;
     
     return 0;
 }
diff --git a/src/libslic3r/MeshBoolean.cpp b/src/libslic3r/MeshBoolean.cpp
index 1c848eb5f..6db1493ec 100644
--- a/src/libslic3r/MeshBoolean.cpp
+++ b/src/libslic3r/MeshBoolean.cpp
@@ -10,17 +10,24 @@
 #include <CGAL/Polygon_mesh_processing/corefinement.h>
 #include <CGAL/Exact_integer.h>
 #include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/repair.h>
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+#include <CGAL/Cartesian_converter.h>
 
 namespace Slic3r {
 namespace MeshBoolean {
 
-typedef Eigen::Map<const Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>> MapMatrixXfUnaligned;
-typedef Eigen::Map<const Eigen::Matrix<int,   Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>> MapMatrixXiUnaligned;
+using MapMatrixXfUnaligned = Eigen::Map<const Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>>;
+using MapMatrixXiUnaligned = Eigen::Map<const Eigen::Matrix<int,   Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor | Eigen::DontAlign>>;
 
-typedef std::pair<Eigen::MatrixXd, Eigen::MatrixXi> EigenMesh;
-
-static TriangleMesh eigen_to_triangle_mesh(const Eigen::MatrixXd& VC, const Eigen::MatrixXi& FC)
+TriangleMesh eigen_to_triangle_mesh(const EigenMesh &emesh)
 {
+    auto &VC = emesh.first; auto &FC = emesh.second;
+    
     Pointf3s points(size_t(VC.rows())); 
     std::vector<Vec3crd> facets(size_t(FC.rows()));
     
@@ -35,7 +42,7 @@ static TriangleMesh eigen_to_triangle_mesh(const Eigen::MatrixXd& VC, const Eige
     return out;
 }
 
-static EigenMesh triangle_mesh_to_eigen_mesh(const TriangleMesh &mesh)
+EigenMesh triangle_mesh_to_eigen(const TriangleMesh &mesh)
 {
     EigenMesh emesh;
     emesh.first = MapMatrixXfUnaligned(mesh.its.vertices.front().data(),
@@ -48,70 +55,116 @@ static EigenMesh triangle_mesh_to_eigen_mesh(const TriangleMesh &mesh)
     return emesh;
 }
 
-void minus(TriangleMesh& A, const TriangleMesh& B)
+void minus(EigenMesh &A, const EigenMesh &B)
 {
-    auto [VA, FA] = triangle_mesh_to_eigen_mesh(A);
-    auto [VB, FB] = triangle_mesh_to_eigen_mesh(B);
-
+    auto &[VA, FA] = A;
+    auto &[VB, FB] = B;
+    
     Eigen::MatrixXd VC;
     Eigen::MatrixXi FC;
     igl::MeshBooleanType boolean_type(igl::MESH_BOOLEAN_TYPE_MINUS);
     igl::copyleft::cgal::mesh_boolean(VA, FA, VB, FB, boolean_type, VC, FC);
-
-    A = eigen_to_triangle_mesh(VC, FC);
+    
+    VA = std::move(VC); FA = std::move(FC);
 }
 
-void self_union(TriangleMesh& mesh)
+void minus(TriangleMesh& A, const TriangleMesh& B)
 {
-    auto [V, F] = triangle_mesh_to_eigen_mesh(mesh);
+    EigenMesh eA = triangle_mesh_to_eigen(A);
+    minus(eA, triangle_mesh_to_eigen(B));
+    A = eigen_to_triangle_mesh(eA);
+}
 
-    Eigen::MatrixXd VC;
-    Eigen::MatrixXi FC;
+void self_union(EigenMesh &A)
+{
+    EigenMesh result;
+    auto &[V, F] = A;
+    auto &[VC, FC] = result;
 
     igl::MeshBooleanType boolean_type(igl::MESH_BOOLEAN_TYPE_UNION);
     igl::copyleft::cgal::mesh_boolean(V, F, Eigen::MatrixXd(), Eigen::MatrixXi(), boolean_type, VC, FC);
     
-    mesh = eigen_to_triangle_mesh(VC, FC);
+    A = std::move(result);
+}
+
+void self_union(TriangleMesh& mesh)
+{
+    auto eM = triangle_mesh_to_eigen(mesh);
+    self_union(eM);
+    mesh = eigen_to_triangle_mesh(eM);
 }
 
 namespace cgal {
 
-namespace CGALProc   = CGAL::Polygon_mesh_processing;
-namespace CGALParams = CGAL::Polygon_mesh_processing::parameters;
+namespace CGALProc    = CGAL::Polygon_mesh_processing;
+namespace CGALParams  = CGAL::Polygon_mesh_processing::parameters;
 
-using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using _CGALMesh = CGAL::Surface_mesh<Kernel::Point_3>;
+using EpecKernel = CGAL::Exact_predicates_exact_constructions_kernel;
+using EpicKernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using _EpicMesh = CGAL::Surface_mesh<EpicKernel::Point_3>;
+using _EpecMesh = CGAL::Surface_mesh<EpecKernel::Point_3>;
 
-struct CGALMesh { _CGALMesh m; };
+struct CGALMesh { _EpicMesh m; };
 
-static void triangle_mesh_to_cgal(const TriangleMesh &M, _CGALMesh &out)
+// /////////////////////////////////////////////////////////////////////////////
+// Converions from and to CGAL mesh
+// /////////////////////////////////////////////////////////////////////////////
+
+template<class _Mesh> void triangle_mesh_to_cgal(const TriangleMesh &M, _Mesh &out)
 {
-    for (const Vec3f &v : M.its.vertices)
-        out.add_vertex(_CGALMesh::Point(v.x(), v.y(), v.z()));
-
-    for (const Vec3crd &face : M.its.indices) {
-        auto f = face.cast<CGAL::SM_Vertex_index>();
-        out.add_face(f(0), f(1), f(2));
+    using Index3 = std::array<size_t, 3>;
+    std::vector<typename _Mesh::Point> points;
+    std::vector<Index3> indices;
+    points.reserve(M.its.vertices.size());
+    indices.reserve(M.its.indices.size());
+    for (auto &v : M.its.vertices) points.emplace_back(v.x(), v.y(), v.z());
+    for (auto &_f : M.its.indices) {
+        auto f = _f.cast<size_t>();
+        indices.emplace_back(Index3{f(0), f(1), f(2)});
     }
+
+    CGALProc::orient_polygon_soup(points, indices);
+    CGALProc::polygon_soup_to_polygon_mesh(points, indices, out);
+    
+    // Number the faces because 'orient_to_bound_a_volume' needs a face <--> index map
+    unsigned index = 0;
+    for (auto face : out.faces()) face = CGAL::SM_Face_index(index++);
+    
+    if(CGAL::is_closed(out))
+        CGALProc::orient_to_bound_a_volume(out);
+    else
+        std::runtime_error("Mesh not watertight");
 }
 
-static TriangleMesh cgal_to_triangle_mesh(const _CGALMesh &cgalmesh)
+inline Vec3d to_vec3d(const _EpicMesh::Point &v)
+{
+    return {v.x(), v.y(), v.z()};
+}
+
+inline Vec3d to_vec3d(const _EpecMesh::Point &v)
+{
+    CGAL::Cartesian_converter<EpecKernel, EpicKernel> cvt;
+    auto iv = cvt(v);
+    return {iv.x(), iv.y(), iv.z()};
+}
+
+template<class _Mesh> TriangleMesh cgal_to_triangle_mesh(const _Mesh &cgalmesh)
 {
     Pointf3s points;
     std::vector<Vec3crd> facets;
     points.reserve(cgalmesh.num_vertices());
     facets.reserve(cgalmesh.num_faces());
-
+    
     for (auto &vi : cgalmesh.vertices()) {
         auto &v = cgalmesh.point(vi); // Don't ask...
-        points.emplace_back(v.x(), v.y(), v.z());
+        points.emplace_back(to_vec3d(v));
     }
-
+    
     for (auto &face : cgalmesh.faces()) {
         auto    vtc = cgalmesh.vertices_around_face(cgalmesh.halfedge(face));
         int     i   = 0;
         Vec3crd trface;
-        for (auto v : vtc) trface(i++) = static_cast<unsigned>(v);
+        for (auto v : vtc) trface(i++) = static_cast<int>(v);
         facets.emplace_back(trface);
     }
     
@@ -120,59 +173,100 @@ static TriangleMesh cgal_to_triangle_mesh(const _CGALMesh &cgalmesh)
     return out;
 }
 
-std::unique_ptr<CGALMesh> triangle_mesh_to_cgal(const TriangleMesh &M)
+std::unique_ptr<CGALMesh, CGALMeshDeleter> triangle_mesh_to_cgal(const TriangleMesh &M)
 {
-    auto out = std::make_unique<CGALMesh>();
+    std::unique_ptr<CGALMesh, CGALMeshDeleter> out(new CGALMesh{});
     triangle_mesh_to_cgal(M, out->m);
     return out;
 }
 
-void cgal_to_triangle_mesh(const CGALMesh &cgalmesh, TriangleMesh &out)
+TriangleMesh cgal_to_triangle_mesh(const CGALMesh &cgalmesh)
 {
-    out = cgal_to_triangle_mesh(cgalmesh.m);
+    return cgal_to_triangle_mesh(cgalmesh.m);
 }
 
-void minus(CGALMesh &A, CGALMesh &B)
+// /////////////////////////////////////////////////////////////////////////////
+// Boolean operations for CGAL meshes
+// /////////////////////////////////////////////////////////////////////////////
+
+static bool _cgal_diff(CGALMesh &A, CGALMesh &B, CGALMesh &R)
 {
-    CGALProc::corefine_and_compute_difference(A.m, B.m, A.m);
+    const auto &p = CGALParams::throw_on_self_intersection(true);
+    return CGALProc::corefine_and_compute_difference(A.m, B.m, R.m, p, p);
 }
 
-void self_union(CGALMesh &A)
+static bool _cgal_union(CGALMesh &A, CGALMesh &B, CGALMesh &R)
 {
-    CGALProc::corefine(A.m, A.m);
+    const auto &p = CGALParams::throw_on_self_intersection(true);
+    return CGALProc::corefine_and_compute_union(A.m, B.m, R.m, p, p);
 }
 
-void minus(TriangleMesh &A, const TriangleMesh &B)
-{   
+static bool _cgal_intersection(CGALMesh &A, CGALMesh &B, CGALMesh &R)
+{
+    const auto &p = CGALParams::throw_on_self_intersection(true);
+    return CGALProc::corefine_and_compute_intersection(A.m, B.m, R.m, p, p);
+}
+
+template<class Op> void _cgal_do(Op &&op, CGALMesh &A, CGALMesh &B)
+{
+    bool success = false;
+    try {
+        CGALMesh result;
+        success = op(A, B, result);
+        A = std::move(result);      // In-place operation does not work
+    } catch (...) {
+        success = false;
+    }
+
+    if (! success)
+        throw std::runtime_error("CGAL mesh boolean operation failed.");
+}
+
+void minus(CGALMesh &A, CGALMesh &B) { _cgal_do(_cgal_diff, A, B); }
+void plus(CGALMesh &A, CGALMesh &B) { _cgal_do(_cgal_union, A, B); }
+void intersect(CGALMesh &A, CGALMesh &B) { _cgal_do(_cgal_intersection, A, B); }
+bool does_self_intersect(const CGALMesh &mesh) { return CGALProc::does_self_intersect(mesh.m); }
+
+// /////////////////////////////////////////////////////////////////////////////
+// Now the public functions for TriangleMesh input:
+// /////////////////////////////////////////////////////////////////////////////
+
+template<class Op> void _mesh_boolean_do(Op &&op, TriangleMesh &A, const TriangleMesh &B)
+{
     CGALMesh meshA;
     CGALMesh meshB;
     triangle_mesh_to_cgal(A, meshA.m);
     triangle_mesh_to_cgal(B, meshB.m);
     
-    CGALMesh meshResult;
-    bool success = false;
-    try {
-        success = CGALProc::corefine_and_compute_difference(meshA.m, meshB.m, meshResult.m,
-            CGALParams::throw_on_self_intersection(true), CGALParams::throw_on_self_intersection(true));
-    }
-    catch (const CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception&) {
-        success = false;
-    }
-    if (! success)
-        throw std::runtime_error("CGAL corefine_and_compute_difference failed");
-
-    A = cgal_to_triangle_mesh(meshResult.m);
-}
-
-void self_union(TriangleMesh &m)
-{
-    _CGALMesh cgalmesh;
-    triangle_mesh_to_cgal(m, cgalmesh);
-    CGALProc::corefine(cgalmesh, cgalmesh);
+    _cgal_do(op, meshA, meshB);
     
-    m = cgal_to_triangle_mesh(cgalmesh);
+    A = cgal_to_triangle_mesh(meshA.m);
 }
 
+void minus(TriangleMesh &A, const TriangleMesh &B)
+{
+    _mesh_boolean_do(_cgal_diff, A, B);
+}
+
+void plus(TriangleMesh &A, const TriangleMesh &B)
+{
+    _mesh_boolean_do(_cgal_union, A, B);
+}
+
+void intersect(TriangleMesh &A, const TriangleMesh &B)
+{
+    _mesh_boolean_do(_cgal_intersection, A, B);
+}
+
+bool does_self_intersect(const TriangleMesh &mesh)
+{
+    CGALMesh cgalm;
+    triangle_mesh_to_cgal(mesh, cgalm.m);
+    return CGALProc::does_self_intersect(cgalm.m);
+}
+
+void CGALMeshDeleter::operator()(CGALMesh *ptr) { delete ptr; }
+
 } // namespace cgal
 
 } // namespace MeshBoolean
diff --git a/src/libslic3r/MeshBoolean.hpp b/src/libslic3r/MeshBoolean.hpp
index 14e3d3b7b..5298906b4 100644
--- a/src/libslic3r/MeshBoolean.hpp
+++ b/src/libslic3r/MeshBoolean.hpp
@@ -2,13 +2,23 @@
 #define libslic3r_MeshBoolean_hpp_
 
 #include <memory>
+#include <exception>
+
+#include <libslic3r/TriangleMesh.hpp>
+#include <Eigen/Geometry>
 
 namespace Slic3r {
 
-class TriangleMesh;
-
 namespace MeshBoolean {
 
+using EigenMesh = std::pair<Eigen::MatrixXd, Eigen::MatrixXi>;
+
+TriangleMesh eigen_to_triangle_mesh(const EigenMesh &emesh);
+EigenMesh triangle_mesh_to_eigen(const TriangleMesh &mesh);
+
+void minus(EigenMesh &A, const EigenMesh &B);
+void self_union(EigenMesh &A);
+    
 void minus(TriangleMesh& A, const TriangleMesh& B);
 void self_union(TriangleMesh& mesh);
 
@@ -16,20 +26,22 @@ namespace cgal {
 
 struct CGALMesh;
 
-std::unique_ptr<CGALMesh> triangle_mesh_to_cgal(const TriangleMesh &M);
-void cgal_to_triangle_mesh(const CGALMesh &cgalmesh, TriangleMesh &out);
+struct CGALMeshDeleter { void operator()(CGALMesh *ptr); };
+
+std::unique_ptr<CGALMesh, CGALMeshDeleter> triangle_mesh_to_cgal(const TriangleMesh &M);
+TriangleMesh cgal_to_triangle_mesh(const CGALMesh &cgalmesh);
     
 // Do boolean mesh difference with CGAL bypassing igl.
 void minus(TriangleMesh &A, const TriangleMesh &B);
+void plus(TriangleMesh &A, const TriangleMesh &B);
+void intersect(TriangleMesh &A, const TriangleMesh &B);
 
-// Do self union only with CGAL.
-void self_union(TriangleMesh& mesh);
-
-// does A = A - B
-// CGAL takes non-const objects as arguments. I suppose it doesn't change B but
-// there is no official garantee.
 void minus(CGALMesh &A, CGALMesh &B);
-void self_union(CGALMesh &A);
+void plus(CGALMesh &A, CGALMesh &B);
+void intersect(CGALMesh &A, CGALMesh &B);
+
+bool does_self_intersect(const TriangleMesh &mesh);
+bool does_self_intersect(const CGALMesh &mesh);
 
 }
 
diff --git a/src/libslic3r/SLAPrintSteps.cpp b/src/libslic3r/SLAPrintSteps.cpp
index 7810a351f..0c293c7fb 100644
--- a/src/libslic3r/SLAPrintSteps.cpp
+++ b/src/libslic3r/SLAPrintSteps.cpp
@@ -143,11 +143,13 @@ void SLAPrint::Steps::drill_holes(SLAPrintObject &po)
         holes_mesh.merge(sla::to_triangle_mesh(holept.to_mesh()));
     
     holes_mesh.require_shared_vertices();
-    MeshBoolean::self_union(holes_mesh);
+    if (!holes_mesh.is_manifold() || MeshBoolean::cgal::does_self_intersect(holes_mesh)) {
+        MeshBoolean::self_union(holes_mesh);
+    }
     
     try {
         MeshBoolean::cgal::minus(hollowed_mesh, holes_mesh);
-    } catch (const std::runtime_error&) {
+    } catch (const std::runtime_error &) {
         throw std::runtime_error(L(
             "Drilling holes into the mesh failed. "
             "This is usually caused by broken model. Try to fix it first."));
diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt
index adcb2722d..97d9b0bc7 100644
--- a/tests/libslic3r/CMakeLists.txt
+++ b/tests/libslic3r/CMakeLists.txt
@@ -12,6 +12,7 @@ add_executable(${_TEST_NAME}_tests
 	test_polygon.cpp
 	test_stl.cpp
         test_meshsimplify.cpp
+        test_meshboolean.cpp
 	)
 
 if (TARGET OpenVDB::openvdb)
diff --git a/tests/libslic3r/test_meshboolean.cpp b/tests/libslic3r/test_meshboolean.cpp
new file mode 100644
index 000000000..97d03ac23
--- /dev/null
+++ b/tests/libslic3r/test_meshboolean.cpp
@@ -0,0 +1,26 @@
+#include <catch2/catch.hpp>
+#include <test_utils.hpp>
+
+#include <libslic3r/TriangleMesh.hpp>
+#include <libslic3r/MeshBoolean.hpp>
+#include <libslic3r/SimplifyMesh.hpp>
+
+using namespace Slic3r;
+
+TEST_CASE("CGAL and TriangleMesh conversions", "[MeshBoolean]") {
+    TriangleMesh sphere = make_sphere(1.);
+    
+    auto cgalmesh_ptr = MeshBoolean::cgal::triangle_mesh_to_cgal(sphere);
+    
+    REQUIRE(cgalmesh_ptr);
+    REQUIRE(! MeshBoolean::cgal::does_self_intersect(*cgalmesh_ptr));
+    
+    TriangleMesh M = MeshBoolean::cgal::cgal_to_triangle_mesh(*cgalmesh_ptr);
+    
+    REQUIRE(M.its.vertices.size() == sphere.its.vertices.size());
+    REQUIRE(M.its.indices.size() == sphere.its.indices.size());
+    
+    REQUIRE(M.volume() == Approx(sphere.volume()));
+    
+    REQUIRE(! MeshBoolean::cgal::does_self_intersect(M));
+}