NEW: add boolean in export stl

Change-Id: Ibeab33f27ad7a2531bb256edd0e3f853cf8def7c
(cherry picked from commit 6cfb4371a1b9053b54569d91a9c29af97a753c0e)

src/libslic3r/CSGMesh/CSGMesh.hpp

@@ -0,0 +1,86 @@
+#ifndef CSGMESH_HPP
+#define CSGMESH_HPP
+
+#include <libslic3r/AnyPtr.hpp>
+#include <admesh/stl.h>
+
+namespace Slic3r { namespace csg {
+
+// A CSGPartT should be an object that can provide at least a mesh + trafo and an
+// associated csg operation. A collection of CSGPartT objects can then
+// be interpreted as one model and used in various contexts. It can be assembled
+// with CGAL or OpenVDB, rendered with OpenCSG or provided to a ray-tracer to
+// deal with various parts of it according to the supported CSG types...
+//
+// A few simple templated interface functions are provided here and a default
+// CSGPart class that implements the necessary means to be usable as a
+// CSGPartT object.
+
+// Supported CSG operation types
+enum class CSGType { Union, Difference, Intersection };
+
+// A CSG part can instruct the processing to push the sub-result until a new
+// csg part with a pop instruction appears. This can be used to implement
+// parentheses in a CSG expression represented by the collection of csg parts.
+// A CSG part can not contain another CSG collection, only a mesh, this is why
+// its easier to do this stacking instead of recursion in the data definition.
+// CSGStackOp::Continue means no stack operation required.
+// When a CSG part contains a Push instruction, it is expected that the CSG
+// operation it contains refers to the whole collection spanning to the nearest
+// part with a Pop instruction.
+// e.g.:
+// {
+//      CUBE1: { mesh: cube, op: Union, stack op: Continue },
+//      CUBE2: { mesh: cube, op: Difference, stack op: Push},
+//      CUBE3: { mesh: cube, op: Union, stack op: Pop}
+// }
+// is a collection of csg parts representing the expression CUBE1 - (CUBE2 + CUBE3)
+enum class CSGStackOp { Push, Continue, Pop };
+
+// Get the CSG operation of the part. Can be overriden for any type
+template<class CSGPartT> CSGType get_operation(const CSGPartT &part)
+{
+    return part.operation;
+}
+
+// Get the stack operation required by the CSG part.
+template<class CSGPartT> CSGStackOp get_stack_operation(const CSGPartT &part)
+{
+    return part.stack_operation;
+}
+
+// Get the mesh for the part. Can be overriden for any type
+template<class CSGPartT>
+const indexed_triangle_set *get_mesh(const CSGPartT &part)
+{
+    return part.its_ptr.get();
+}
+
+// Get the transformation associated with the mesh inside a CSGPartT object.
+// Can be overriden for any type.
+template<class CSGPartT>
+Transform3f get_transform(const CSGPartT &part)
+{
+    return part.trafo;
+}
+
+// Default implementation
+struct CSGPart {
+    AnyPtr<const indexed_triangle_set> its_ptr;
+    Transform3f trafo;
+    CSGType operation;
+    CSGStackOp stack_operation;
+
+    CSGPart(AnyPtr<const indexed_triangle_set> ptr = {},
+            CSGType                            op  = CSGType::Union,
+            const Transform3f                 &tr  = Transform3f::Identity())
+        : its_ptr{std::move(ptr)}
+        , operation{op}
+        , stack_operation{CSGStackOp::Continue}
+        , trafo{tr}
+    {}
+};
+
+}} // namespace Slic3r::csg
+
+#endif // CSGMESH_HPP

src/libslic3r/CSGMesh/CSGMeshCopy.hpp

@@ -0,0 +1,80 @@
+#ifndef CSGMESHCOPY_HPP
+#define CSGMESHCOPY_HPP
+
+#include "CSGMesh.hpp"
+
+namespace Slic3r { namespace csg {
+
+// Copy a csg range but for the meshes, only copy the pointers. If the copy
+// is made from a CSGPart compatible object, and the pointer is a shared one,
+// it will be copied with reference counting.
+template<class It, class OutIt>
+void copy_csgrange_shallow(const Range<It> &csgrange, OutIt out)
+{
+    for (const auto &part : csgrange) {
+        CSGPart cpy{{},
+                    get_operation(part),
+                    get_transform(part)};
+
+        cpy.stack_operation = get_stack_operation(part);
+
+        if constexpr (std::is_convertible_v<decltype(part), const CSGPart&>) {
+            if (auto shptr = part.its_ptr.get_shared_cpy()) {
+                cpy.its_ptr = shptr;
+            }
+        }
+
+        if (!cpy.its_ptr)
+            cpy.its_ptr = AnyPtr<const indexed_triangle_set>{get_mesh(part)};
+
+        *out = std::move(cpy);
+        ++out;
+    }
+}
+
+// Copy the csg range, allocating new meshes
+template<class It, class OutIt>
+void copy_csgrange_deep(const Range<It> &csgrange, OutIt out)
+{
+    for (const auto &part : csgrange) {
+
+        CSGPart cpy{{}, get_operation(part), get_transform(part)};
+
+        if (auto meshptr = get_mesh(part)) {
+            cpy.its_ptr = std::make_unique<const indexed_triangle_set>(*meshptr);
+        }
+
+        cpy.stack_operation = get_stack_operation(part);
+
+        *out = std::move(cpy);
+        ++out;
+    }
+}
+
+template<class ItA, class ItB>
+bool is_same(const Range<ItA> &A, const Range<ItB> &B)
+{
+    bool ret = true;
+
+    size_t s = A.size();
+
+    if (B.size() != s)
+        ret = false;
+
+    size_t i = 0;
+    auto itA = A.begin();
+    auto itB = B.begin();
+    for (; ret && i < s; ++itA, ++itB, ++i) {
+        ret = ret &&
+              get_mesh(*itA) == get_mesh(*itB) &&
+              get_operation(*itA) == get_operation(*itB) &&
+              get_stack_operation(*itA) == get_stack_operation(*itB) &&
+              get_transform(*itA).isApprox(get_transform(*itB));
+    }
+
+    return ret;
+}
+
+}} // namespace Slic3r::csg
+
+#endif // CSGCOPY_HPP

src/libslic3r/CSGMesh/ModelToCSGMesh.hpp

@@ -0,0 +1,88 @@
+#ifndef MODELTOCSGMESH_HPP
+#define MODELTOCSGMESH_HPP
+
+#include "CSGMesh.hpp"
+
+#include "libslic3r/Model.hpp"
+#include "libslic3r/SLA/Hollowing.hpp"
+#include "libslic3r/MeshSplitImpl.hpp"
+
+namespace Slic3r { namespace csg {
+
+// Flags to select which parts to export from Model into a csg part collection.
+// These flags can be chained with the | operator
+enum ModelParts {
+    mpartsPositive = 1,   // Include positive parts
+    mpartsNegative = 2,   // Include negative parts
+    mpartsDrillHoles = 4, // Include drill holes
+    mpartsDoSplits = 8,   // Split each splitable mesh and export as a union of csg parts
+};
+
+template<class OutIt>
+void model_to_csgmesh(const ModelObject &mo,
+                      const Transform3d &trafo, // Applies to all exported parts
+                      OutIt              out,   // Output iterator
+                      // values of ModelParts OR-ed
+                      int                parts_to_include = mpartsPositive
+                      )
+{
+    bool do_positives  = parts_to_include & mpartsPositive;
+    bool do_negatives  = parts_to_include & mpartsNegative;
+    bool do_drillholes = parts_to_include & mpartsDrillHoles;
+    bool do_splits     = parts_to_include & mpartsDoSplits;
+
+    for (const ModelVolume *vol : mo.volumes) {
+        if (vol && vol->mesh_ptr() &&
+            ((do_positives && vol->is_model_part()) ||
+             (do_negatives && vol->is_negative_volume()))) {
+
+            if (do_splits && its_is_splittable(vol->mesh().its)) {
+                CSGPart part_begin{{}, vol->is_model_part() ? CSGType::Union : CSGType::Difference};
+                part_begin.stack_operation = CSGStackOp::Push;
+                *out = std::move(part_begin);
+                ++out;
+
+                its_split(vol->mesh().its, SplitOutputFn{[&out, &vol, &trafo](indexed_triangle_set &&its) {
+                              if (its.empty())
+                                  return;
+
+                              CSGPart part{std::make_unique<indexed_triangle_set>(std::move(its)),
+                                       CSGType::Union,
+                                       (trafo * vol->get_matrix()).cast<float>()};
+
+                              *out = std::move(part);
+                              ++out;
+                          }});
+
+                CSGPart part_end{{}};
+                part_end.stack_operation = CSGStackOp::Pop;
+                *out = std::move(part_end);
+                ++out;
+            } else {
+                CSGPart part{&(vol->mesh().its),
+                             vol->is_model_part() ? CSGType::Union : CSGType::Difference,
+                             (trafo * vol->get_matrix()).cast<float>()};
+
+                *out = std::move(part);
+                ++out;
+            }
+        }
+    }
+
+    //if (do_drillholes) {
+    //    sla::DrainHoles drainholes = sla::transformed_drainhole_points(mo, trafo);
+
+    //    for (const sla::DrainHole &dhole : drainholes) {
+    //        CSGPart part{std::make_unique<const indexed_triangle_set>(
+    //                         dhole.to_mesh()),
+    //                     CSGType::Difference};
+
+    //        *out = std::move(part);
+    //        ++out;
+    //    }
+    //}
+}
+
+}} // namespace Slic3r::csg
+
+#endif // MODELTOCSGMESH_HPP

src/libslic3r/CSGMesh/PerformCSGMeshBooleans.hpp

@@ -0,0 +1,205 @@
+#ifndef PERFORMCSGMESHBOOLEANS_HPP
+#define PERFORMCSGMESHBOOLEANS_HPP
+
+#include <stack>
+#include <vector>
+
+#include "CSGMesh.hpp"
+
+#include "libslic3r/Execution/ExecutionTBB.hpp"
+//#include "libslic3r/Execution/ExecutionSeq.hpp"
+#include "libslic3r/MeshBoolean.hpp"
+
+namespace Slic3r { namespace csg {
+
+// This method can be overriden when a specific CSGPart type supports caching
+// of the voxel grid
+template<class CSGPartT>
+MeshBoolean::cgal::CGALMeshPtr get_cgalmesh(const CSGPartT &csgpart)
+{
+    const indexed_triangle_set *its = csg::get_mesh(csgpart);
+    indexed_triangle_set dummy;
+
+    if (!its)
+        its = &dummy;
+
+    MeshBoolean::cgal::CGALMeshPtr ret;
+
+    indexed_triangle_set m = *its;
+    its_transform(m, get_transform(csgpart), true);
+
+    try {
+        ret = MeshBoolean::cgal::triangle_mesh_to_cgal(m);
+    } catch (...) {
+        // errors are ignored, simply return null
+        ret = nullptr;
+    }
+
+    return ret;
+}
+
+
+namespace detail_cgal {
+
+using MeshBoolean::cgal::CGALMeshPtr;
+
+inline void perform_csg(CSGType op, CGALMeshPtr &dst, CGALMeshPtr &src)
+{
+    if (!dst && op == CSGType::Union && src) {
+        dst = std::move(src);
+        return;
+    }
+
+    if (!dst || !src)
+        return;
+
+    switch (op) {
+    case CSGType::Union:
+        MeshBoolean::cgal::plus(*dst, *src);
+        break;
+    case CSGType::Difference:
+        MeshBoolean::cgal::minus(*dst, *src);
+        break;
+    case CSGType::Intersection:
+        MeshBoolean::cgal::intersect(*dst, *src);
+        break;
+    }
+}
+
+template<class Ex, class It>
+std::vector<CGALMeshPtr> get_cgalptrs(Ex policy, const Range<It> &csgrange)
+{
+    std::vector<CGALMeshPtr> ret(csgrange.size());
+    execution::for_each(policy, size_t(0), csgrange.size(),
+                        [&csgrange, &ret](size_t i) {
+        auto it = csgrange.begin();
+        std::advance(it, i);
+        auto &csgpart = *it;
+        ret[i]        = get_cgalmesh(csgpart);
+    });
+
+    return ret;
+}
+
+} // namespace detail
+
+// Process the sequence of CSG parts with CGAL.
+template<class It>
+void perform_csgmesh_booleans_cgal(MeshBoolean::cgal::CGALMeshPtr &cgalm,
+                              const Range<It>                &csgrange)
+{
+    using MeshBoolean::cgal::CGALMesh;
+    using MeshBoolean::cgal::CGALMeshPtr;
+    using namespace detail_cgal;
+
+    struct Frame {
+        CSGType op; CGALMeshPtr cgalptr;
+        explicit Frame(CSGType csgop = CSGType::Union)
+            : op{ csgop }
+            , cgalptr{ MeshBoolean::cgal::triangle_mesh_to_cgal(indexed_triangle_set{}) }
+        {}
+    };
+
+    std::stack opstack{ std::vector<Frame>{} };
+
+    opstack.push(Frame{});
+
+    std::vector<CGALMeshPtr> cgalmeshes = get_cgalptrs(ex_tbb, csgrange);
+
+    size_t csgidx = 0;
+    for (auto& csgpart : csgrange) {
+
+        auto op = get_operation(csgpart);
+        CGALMeshPtr& cgalptr = cgalmeshes[csgidx++];
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Push) {
+            opstack.push(Frame{ op });
+            op = CSGType::Union;
+        }
+
+        Frame* top = &opstack.top();
+
+        perform_csg(get_operation(csgpart), top->cgalptr, cgalptr);
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Pop) {
+            CGALMeshPtr src = std::move(top->cgalptr);
+            auto popop = opstack.top().op;
+            opstack.pop();
+            CGALMeshPtr& dst = opstack.top().cgalptr;
+            perform_csg(popop, dst, src);
+        }
+    }
+
+    cgalm = std::move(opstack.top().cgalptr);
+}
+
+template<class It, class Visitor>
+It check_csgmesh_booleans(const Range<It> &csgrange, Visitor &&vfn)
+{
+    using namespace detail_cgal;
+
+    std::vector<CGALMeshPtr> cgalmeshes(csgrange.size());
+    auto check_part = [&csgrange, &cgalmeshes](size_t i)
+    {
+        auto it = csgrange.begin();
+        std::advance(it, i);
+        auto &csgpart = *it;
+        auto m = get_cgalmesh(csgpart);
+
+        // mesh can be nullptr if this is a stack push or pull
+        if (!get_mesh(csgpart) && get_stack_operation(csgpart) != CSGStackOp::Continue) {
+            cgalmeshes[i] = MeshBoolean::cgal::triangle_mesh_to_cgal(indexed_triangle_set{});
+            return;
+        }
+
+        try {
+            if (!m || MeshBoolean::cgal::empty(*m))
+                return;
+
+            if (!MeshBoolean::cgal::does_bound_a_volume(*m))
+                return;
+
+            if (MeshBoolean::cgal::does_self_intersect(*m))
+                return;
+        }
+        catch (...) { return; }
+
+        cgalmeshes[i] = std::move(m);
+    };
+    execution::for_each(ex_tbb, size_t(0), csgrange.size(), check_part);
+
+    It ret = csgrange.end();
+    for (size_t i = 0; i < csgrange.size(); ++i) {
+        if (!cgalmeshes[i]) {
+            auto it = csgrange.begin();
+            std::advance(it, i);
+            vfn(it);
+
+            if (ret == csgrange.end())
+                ret = it;
+        }
+    }
+
+    return ret;
+}
+
+template<class It>
+It check_csgmesh_booleans(const Range<It> &csgrange)
+{
+        return check_csgmesh_booleans(csgrange, [](auto &) {});
+}
+
+template<class It>
+MeshBoolean::cgal::CGALMeshPtr perform_csgmesh_booleans(const Range<It> &csgparts)
+{
+    auto ret = MeshBoolean::cgal::triangle_mesh_to_cgal(indexed_triangle_set{});
+    if (ret) {
+        perform_csgmesh_booleans_cgal(ret, csgparts);
+    }
+    return ret;
+}
+
+} // namespace csg
+} // namespace Slic3r
+
+#endif // PERFORMCSGMESHBOOLEANS_HPP

src/libslic3r/CSGMesh/SliceCSGMesh.hpp

@@ -0,0 +1,131 @@
+#ifndef SLICECSGMESH_HPP
+#define SLICECSGMESH_HPP
+
+#include "CSGMesh.hpp"
+
+#include <stack>
+
+#include "libslic3r/TriangleMeshSlicer.hpp"
+#include "libslic3r/ClipperUtils.hpp"
+#include "libslic3r/Execution/ExecutionTBB.hpp"
+
+namespace Slic3r { namespace csg {
+
+namespace detail {
+
+inline void merge_slices(csg::CSGType op, size_t i,
+                  std::vector<ExPolygons> &target,
+                  std::vector<ExPolygons> &source)
+{
+    switch(op) {
+    case CSGType::Union:
+        for (ExPolygon &expoly : source[i])
+            target[i].emplace_back(std::move(expoly));
+        break;
+    case CSGType::Difference:
+        target[i] = diff_ex(target[i], source[i]);
+        break;
+    case CSGType::Intersection:
+        target[i] = intersection_ex(target[i], source[i]);
+        break;
+    }
+}
+
+inline void collect_nonempty_indices(csg::CSGType                   op,
+                              const std::vector<float>      &slicegrid,
+                              const std::vector<ExPolygons> &slices,
+                              std::vector<size_t>           &indices)
+{
+    indices.clear();
+    for (size_t i = 0; i < slicegrid.size(); ++i) {
+        if (op == CSGType::Intersection || !slices[i].empty())
+            indices.emplace_back(i);
+    }
+}
+
+} // namespace detail
+
+template<class ItCSG>
+std::vector<ExPolygons> slice_csgmesh_ex(
+    const Range<ItCSG>          &csgrange,
+    const std::vector<float>    &slicegrid,
+    const MeshSlicingParamsEx   &params,
+    const std::function<void()> &throw_on_cancel = [] {})
+{
+    using namespace detail;
+
+    struct Frame { CSGType op; std::vector<ExPolygons> slices; };
+
+    std::stack opstack{std::vector<Frame>{}};
+
+    MeshSlicingParamsEx params_cpy = params;
+    auto trafo = params.trafo;
+    auto nonempty_indices = reserve_vector<size_t>(slicegrid.size());
+
+    opstack.push({CSGType::Union, std::vector<ExPolygons>(slicegrid.size())});
+
+    for (const auto &csgpart : csgrange) {
+        const indexed_triangle_set *its = csg::get_mesh(csgpart);
+
+        auto op = get_operation(csgpart);
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Push) {
+            opstack.push({op, std::vector<ExPolygons>(slicegrid.size())});
+            op = CSGType::Union;
+        }
+
+        Frame *top = &opstack.top();
+
+        if (its) {
+            params_cpy.trafo = trafo * csg::get_transform(csgpart).template cast<double>();
+            std::vector<ExPolygons> slices = slice_mesh_ex(*its,
+                                                           slicegrid, params_cpy,
+                                                           throw_on_cancel);
+
+            assert(slices.size() == slicegrid.size());
+
+            collect_nonempty_indices(op, slicegrid, slices, nonempty_indices);
+
+            execution::for_each(
+                ex_tbb, nonempty_indices.begin(), nonempty_indices.end(),
+                [op, &slices, &top](size_t i) {
+                    merge_slices(op, i, top->slices, slices);
+                }, execution::max_concurrency(ex_tbb));
+        }
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Pop) {
+            std::vector<ExPolygons> popslices = std::move(top->slices);
+            auto popop = opstack.top().op;
+            opstack.pop();
+            std::vector<ExPolygons> &prev_slices = opstack.top().slices;
+
+            collect_nonempty_indices(popop, slicegrid, popslices, nonempty_indices);
+
+            execution::for_each(
+                ex_tbb, nonempty_indices.begin(), nonempty_indices.end(),
+                [&popslices, &prev_slices, popop](size_t i) {
+                    merge_slices(popop, i, prev_slices, popslices);
+                }, execution::max_concurrency(ex_tbb));
+        }
+    }
+
+    std::vector<ExPolygons> ret = std::move(opstack.top().slices);
+
+    // TODO: verify if this part can be omitted or not.
+    execution::for_each(ex_tbb, ret.begin(), ret.end(), [](ExPolygons &slice) {
+        auto it = std::remove_if(slice.begin(), slice.end(), [](const ExPolygon &p){
+            return p.area() < double(SCALED_EPSILON) * double(SCALED_EPSILON);
+        });
+
+        // Hopefully, ExPolygons are moved, not copied to new positions
+        // and that is cheap for expolygons
+        slice.erase(it, slice.end());
+        slice = union_ex(slice);
+    }, execution::max_concurrency(ex_tbb));
+
+    return ret;
+}
+
+}} // namespace Slic3r::csg
+
+#endif // SLICECSGMESH_HPP

src/libslic3r/CSGMesh/TriangleMeshAdapter.hpp

@@ -0,0 +1,95 @@
+#ifndef TRIANGLEMESHADAPTER_HPP
+#define TRIANGLEMESHADAPTER_HPP
+
+#include "CSGMesh.hpp"
+
+#include "libslic3r/TriangleMesh.hpp"
+
+namespace Slic3r { namespace csg {
+
+// Provide default overloads for indexed_triangle_set to be usable as a plain
+// CSGPart with an implicit union operation
+
+inline CSGType get_operation(const indexed_triangle_set &part)
+{
+    return CSGType::Union;
+}
+
+inline CSGStackOp get_stack_operation(const indexed_triangle_set &part)
+{
+    return CSGStackOp::Continue;
+}
+
+inline const indexed_triangle_set * get_mesh(const indexed_triangle_set &part)
+{
+    return ∂
+}
+
+inline Transform3f get_transform(const indexed_triangle_set &part)
+{
+    return Transform3f::Identity();
+}
+
+inline CSGType get_operation(const indexed_triangle_set *const part)
+{
+    return CSGType::Union;
+}
+
+inline CSGStackOp get_stack_operation(const indexed_triangle_set *const part)
+{
+    return CSGStackOp::Continue;
+}
+
+inline const indexed_triangle_set * get_mesh(const indexed_triangle_set *const part)
+{
+    return part;
+}
+
+inline Transform3f get_transform(const indexed_triangle_set *const part)
+{
+    return Transform3f::Identity();
+}
+
+inline CSGType get_operation(const TriangleMesh &part)
+{
+    return CSGType::Union;
+}
+
+inline CSGStackOp get_stack_operation(const TriangleMesh &part)
+{
+    return CSGStackOp::Continue;
+}
+
+inline const indexed_triangle_set * get_mesh(const TriangleMesh &part)
+{
+    return &part.its;
+}
+
+inline Transform3f get_transform(const TriangleMesh &part)
+{
+    return Transform3f::Identity();
+}
+
+inline CSGType get_operation(const TriangleMesh * const part)
+{
+    return CSGType::Union;
+}
+
+inline CSGStackOp get_stack_operation(const TriangleMesh * const part)
+{
+    return CSGStackOp::Continue;
+}
+
+inline const indexed_triangle_set * get_mesh(const TriangleMesh * const part)
+{
+    return &part->its;
+}
+
+inline Transform3f get_transform(const TriangleMesh * const part)
+{
+    return Transform3f::Identity();
+}
+
+}} // namespace Slic3r::csg
+
+#endif // TRIANGLEMESHADAPTER_HPP

src/libslic3r/CSGMesh/VoxelizeCSGMesh.hpp

@@ -0,0 +1,116 @@
+#ifndef VOXELIZECSGMESH_HPP
+#define VOXELIZECSGMESH_HPP
+
+#include <functional>
+#include <stack>
+
+#include "CSGMesh.hpp"
+#include "libslic3r/OpenVDBUtils.hpp"
+#include "libslic3r/Execution/ExecutionTBB.hpp"
+
+namespace Slic3r { namespace csg {
+
+using VoxelizeParams = MeshToGridParams;
+
+// This method can be overriden when a specific CSGPart type supports caching
+// of the voxel grid
+template<class CSGPartT>
+VoxelGridPtr get_voxelgrid(const CSGPartT &csgpart, VoxelizeParams params)
+{
+    const indexed_triangle_set *its = csg::get_mesh(csgpart);
+    VoxelGridPtr ret;
+
+    params.trafo(params.trafo() * csg::get_transform(csgpart));
+
+    if (its)
+        ret = mesh_to_grid(*its, params);
+
+    return ret;
+}
+
+namespace detail {
+
+inline void perform_csg(CSGType op, VoxelGridPtr &dst, VoxelGridPtr &src)
+{
+    if (!dst || !src)
+        return;
+
+    switch (op) {
+    case CSGType::Union:
+        if (is_grid_empty(*dst) && !is_grid_empty(*src))
+            dst = clone(*src);
+        else
+            grid_union(*dst, *src);
+
+        break;
+    case CSGType::Difference:
+        grid_difference(*dst, *src);
+        break;
+    case CSGType::Intersection:
+        grid_intersection(*dst, *src);
+        break;
+    }
+}
+
+} // namespace detail
+
+template<class It>
+VoxelGridPtr voxelize_csgmesh(const Range<It>      &csgrange,
+                              const VoxelizeParams &params = {})
+{
+    using namespace detail;
+
+    VoxelGridPtr ret;
+
+    std::vector<VoxelGridPtr> grids (csgrange.size());
+
+    execution::for_each(ex_tbb, size_t(0), csgrange.size(), [&](size_t csgidx) {
+        if (params.statusfn() && params.statusfn()(-1))
+            return;
+
+        auto it = csgrange.begin();
+        std::advance(it, csgidx);
+        auto &csgpart = *it;
+        grids[csgidx] = get_voxelgrid(csgpart, params);
+    }, execution::max_concurrency(ex_tbb));
+
+    size_t csgidx = 0;
+    struct Frame { CSGType op = CSGType::Union; VoxelGridPtr grid; };
+    std::stack opstack{std::vector<Frame>{}};
+
+    opstack.push({CSGType::Union, mesh_to_grid({}, params)});
+
+    for (auto &csgpart : csgrange) {
+        if (params.statusfn() && params.statusfn()(-1))
+            break;
+
+        auto &partgrid = grids[csgidx++];
+
+        auto op = get_operation(csgpart);
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Push) {
+            opstack.push({op, mesh_to_grid({}, params)});
+            op = CSGType::Union;
+        }
+
+        Frame *top = &opstack.top();
+
+        perform_csg(get_operation(csgpart), top->grid, partgrid);
+
+        if (get_stack_operation(csgpart) == CSGStackOp::Pop) {
+            VoxelGridPtr popgrid = std::move(top->grid);
+            auto popop = opstack.top().op;
+            opstack.pop();
+            VoxelGridPtr &grid = opstack.top().grid;
+            perform_csg(popop, grid, popgrid);
+        }
+    }
+
+    ret = std::move(opstack.top().grid);
+
+    return ret;
+}
+
+}} // namespace Slic3r::csg
+
+#endif // VOXELIZECSGMESH_HPP