Tile the booleans on layers of many pieces

ClipperLib slows down with the number of edges on a scan line, and a
layer cut through a fine relief has tens of thousands of pieces.
detect_surfaces_type ~50 s at 0.1 mm / 2000k, was ~145.

src/libslic3r/ClipperUtils.cpp

@@ -2,6 +2,8 @@
 #include <numeric>
 #include <unordered_map>
 
+#include <tbb/parallel_for.h>
+
 #include "ClipperUtils.hpp"
 #include "Geometry.hpp"
 #include "ShortestPath.hpp"
@@ -813,6 +815,63 @@ Slic3r::ExPolygons intersection_ex(const Slic3r::Surfaces &subject, const Slic3r
     { return _clipper_ex(ClipperLib::ctIntersection, ClipperUtils::SurfacesProvider(subject), ClipperUtils::SurfacesProvider(clip), do_safety_offset); }
 Slic3r::ExPolygons intersection_ex(const Slic3r::SurfacesPtr &subject, const Slic3r::ExPolygons &clip, ApplySafetyOffset do_safety_offset)
     { return _clipper_ex(ClipperLib::ctIntersection, ClipperUtils::SurfacesPtrProvider(subject), ClipperUtils::ExPolygonsProvider(clip), do_safety_offset); }
+
+static Slic3r::ExPolygons clipper_ex_by_piece(ClipperLib::ClipType clipType, const Slic3r::ExPolygons &subject, const Slic3r::Polygons &clip, ApplySafetyOffset do_safety_offset)
+{
+    // The subject ExPolygons are split into tiles by the centres of their boxes, a few dozen per tile, and each tile is one
+    // ClipperLib call with the clip cut to the box of the tile's ExPolygons.
+    BoundingBox extent;
+    std::vector<BoundingBox> bboxes;
+    bboxes.reserve(subject.size());
+    for (const ExPolygon &expoly : subject) {
+        bboxes.emplace_back(get_extents(expoly));
+        extent.merge(bboxes.back());
+    }
+    const int tiles = std::clamp(int(std::sqrt(double(subject.size()) / 32.)), 1, 32);
+    std::vector<std::vector<size_t>> members(size_t(tiles * tiles));
+    std::vector<BoundingBox>         tile_bboxes(members.size());
+    if (extent.defined) {
+        const Point   size   = extent.size();
+        const coord_t tile_w = std::max<coord_t>(1, size.x() / tiles + 1), tile_h = std::max<coord_t>(1, size.y() / tiles + 1);
+        for (size_t i = 0; i < subject.size(); ++i) {
+            const Point  c    = bboxes[i].center();
+            const size_t tile = size_t(std::clamp(int((c.y() - extent.min.y()) / tile_h), 0, tiles - 1) * tiles +
+                                       std::clamp(int((c.x() - extent.min.x()) / tile_w), 0, tiles - 1));
+            members[tile].emplace_back(i);
+            tile_bboxes[tile].merge(bboxes[i]);
+        }
+    }
+    std::vector<BoundingBox> clip_bboxes;
+    clip_bboxes.reserve(clip.size());
+    for (const Polygon &polygon : clip)
+        clip_bboxes.emplace_back(get_extents(polygon));
+
+    std::vector<Slic3r::ExPolygons> out_tiles(members.size());
+    tbb::parallel_for(size_t(0), members.size(), [&](size_t tile) {
+        if (members[tile].empty())
+            return;
+        Slic3r::ExPolygons local_subject;
+        local_subject.reserve(members[tile].size());
+        for (size_t i : members[tile])
+            local_subject.emplace_back(subject[i]);
+        // Grown so that the cut edges of the clip stay clear of the subject, also after the safety offset.
+        const BoundingBox bbox = tile_bboxes[tile].inflated(SCALED_EPSILON);
+        Polygons local_clip;
+        for (size_t i = 0; i < clip.size(); ++i)
+            if (clip_bboxes[i].overlap(bbox))
+                if (Polygon clipped = ClipperUtils::clip_clipper_polygon_with_subject_bbox(clip[i], bbox); ! clipped.empty())
+                    local_clip.emplace_back(std::move(clipped));
+        out_tiles[tile] = _clipper_ex(clipType, ClipperUtils::ExPolygonsProvider(local_subject), ClipperUtils::PolygonsProvider(local_clip), do_safety_offset);
+    });
+    Slic3r::ExPolygons out;
+    for (Slic3r::ExPolygons &out_tile : out_tiles)
+        append(out, std::move(out_tile));
+    return out;
+}
+Slic3r::ExPolygons diff_ex_by_piece(const Slic3r::ExPolygons &subject, const Slic3r::Polygons &clip, ApplySafetyOffset do_safety_offset)
+    { return clipper_ex_by_piece(ClipperLib::ctDifference, subject, clip, do_safety_offset); }
+Slic3r::ExPolygons intersection_ex_by_piece(const Slic3r::ExPolygons &subject, const Slic3r::Polygons &clip, ApplySafetyOffset do_safety_offset)
+    { return clipper_ex_by_piece(ClipperLib::ctIntersection, subject, clip, do_safety_offset); }
 // May be used to "heal" unusual models (3DLabPrints etc.) by providing fill_type (pftEvenOdd, pftNonZero, pftPositive, pftNegative).
 Slic3r::ExPolygons union_ex(const Slic3r::Polygons &subject, ClipperLib::PolyFillType fill_type)
     { return _clipper_ex(ClipperLib::ctUnion, ClipperUtils::PolygonsProvider(subject), ClipperUtils::EmptyPathsProvider(), ApplySafetyOffset::No, fill_type); }

src/libslic3r/ClipperUtils.hpp

@@ -518,6 +518,11 @@ Slic3r::ExPolygons intersection_ex(const Slic3r::Surfaces &subject, const Slic3r
 Slic3r::ExPolygons intersection_ex(const Slic3r::Surfaces &subject, const Slic3r::ExPolygons &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
 Slic3r::ExPolygons intersection_ex(const Slic3r::Surfaces &subject, const Slic3r::Surfaces &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
 Slic3r::ExPolygons intersection_ex(const Slic3r::SurfacesPtr &subject, const Slic3r::ExPolygons &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
+// diff_ex() / intersection_ex() of the subject split into tiles, each against only the part of the clip near it, the tiles in
+// parallel. The same area as the operation on the whole subject when its ExPolygons do not overlap, and much faster for a
+// subject of thousands of pieces spread over a layer: ClipperLib slows down with the number of edges crossing a scan line.
+Slic3r::ExPolygons diff_ex_by_piece(const Slic3r::ExPolygons &subject, const Slic3r::Polygons &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
+Slic3r::ExPolygons intersection_ex_by_piece(const Slic3r::ExPolygons &subject, const Slic3r::Polygons &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
 Slic3r::Polylines  intersection_pl(const Slic3r::Polylines &subject, const Slic3r::Polygon &clip);
 Slic3r::Polylines  intersection_pl(const Slic3r::Polyline &subject, const Slic3r::ExPolygon &clip);
 Slic3r::Polylines  intersection_pl(const Slic3r::Polylines &subject, const Slic3r::ExPolygon &clip);

src/libslic3r/LayerRegion.cpp

@@ -72,10 +72,11 @@ void LayerRegion::slices_to_fill_surfaces_clipped()
         by_surface[size_t(surface.surface_type)].emplace_back(&surface);
     // Trim surfaces by the fill_boundaries.
     this->fill_surfaces.surfaces.clear();
+    const Polygons fill_boundaries = to_polygons(this->fill_expolygons);
     for (size_t surface_type = 0; surface_type < size_t(stCount); ++ surface_type) {
         const SurfacesPtr &this_surfaces = by_surface[surface_type];
         if (! this_surfaces.empty())
-            this->fill_surfaces.append(intersection_ex(this_surfaces, this->fill_expolygons), SurfaceType(surface_type));
+            this->fill_surfaces.append(intersection_ex_by_piece(to_expolygons(this_surfaces), fill_boundaries), SurfaceType(surface_type));
     }
 }
 

src/libslic3r/MultiMaterialSegmentation.cpp

@@ -1853,53 +1853,6 @@ static void remove_multiple_edges_in_vertices(MMU_Graph &graph, const std::vecto
     }
 }
 
-// diff_ex(subject, clip), one subject ExPolygon at a time against the part of `clip` in its box. The pieces are disjoint,
-// so the result is the same; but ClipperLib never gets a whole finely painted layer at once, where re-linking the holes
-// of its PolyTree (FixupFirstLefts) is quadratic in the thousands of regions.
-static ExPolygons diff_ex_by_piece(const ExPolygons &subject, const ExPolygons &clip)
-{
-    // A coarse grid over the clip's boxes, so each piece only looks at the clip regions near it.
-    std::vector<BoundingBox> clip_bboxes;
-    clip_bboxes.reserve(clip.size());
-    BoundingBox extent;
-    for (const ExPolygon &expoly : clip) {
-        clip_bboxes.emplace_back(get_extents(expoly));
-        extent.merge(clip_bboxes.back());
-    }
-    constexpr int GRID = 64;
-    const Point   size = extent.defined ? extent.size() : Point(1, 1);
-    const coord_t cell_w = std::max<coord_t>(1, size.x() / GRID + 1), cell_h = std::max<coord_t>(1, size.y() / GRID + 1);
-    const auto    cells  = [&](const BoundingBox &bb, auto &&fn) {
-        const int x0 = std::clamp(int((bb.min.x() - extent.min.x()) / cell_w), 0, GRID - 1), x1 = std::clamp(int((bb.max.x() - extent.min.x()) / cell_w), 0, GRID - 1);
-        const int y0 = std::clamp(int((bb.min.y() - extent.min.y()) / cell_h), 0, GRID - 1), y1 = std::clamp(int((bb.max.y() - extent.min.y()) / cell_h), 0, GRID - 1);
-        for (int y = y0; y <= y1; ++y)
-            for (int x = x0; x <= x1; ++x)
-                fn(y * GRID + x);
-    };
-    std::vector<std::vector<size_t>> grid(GRID * GRID);
-    if (extent.defined)
-        for (size_t i = 0; i < clip.size(); ++i)
-            cells(clip_bboxes[i], [&](int cell) { grid[cell].emplace_back(i); });
-
-    std::vector<ExPolygons> pieces(subject.size());
-    tbb::parallel_for(size_t(0), subject.size(), [&](size_t idx) {
-        const BoundingBox bbox = get_extents(subject[idx]).inflated(SCALED_EPSILON);
-        std::vector<size_t> near;
-        if (extent.defined && bbox.overlap(extent)) {
-            cells(bbox, [&](int cell) { append(near, grid[cell]); });
-            sort_remove_duplicates(near);
-        }
-        Polygons nearby_clip;
-        for (size_t i : near)
-            if (clip_bboxes[i].overlap(bbox))
-                polygons_append(nearby_clip, ClipperUtils::clip_clipper_polygons_with_subject_bbox(clip[i], bbox));
-        pieces[idx] = nearby_clip.empty() ? ExPolygons{ subject[idx] } : diff_ex(subject[idx], nearby_clip);
-    });
-    ExPolygons out;
-    for (ExPolygons &piece : pieces)
-        append(out, std::move(piece));
-    return out;
-}
 
 // Finds the islands (layer ExPolygons) a region piece overlaps. A top or bottom region is projected from the neighbouring
 // layers and may reach past the island it belongs to, or over several islands.
@@ -2031,9 +1984,9 @@ static std::vector<std::vector<ExPolygons>> merge_segmented_layers(const std::ve
             // Side regions minus the top/bottom regions of every colour.
             std::vector<std::vector<ExPolygons>> merged(num_buckets, std::vector<ExPolygons>(num_facets_states));
             tbb::parallel_for(size_t(0), num_buckets, [&](size_t bucket) {
-                ExPolygons tops_all;
+                Polygons tops_all;
                 for (const ExPolygons &t : tops[bucket])
-                    append(tops_all, t);
+                    polygons_append(tops_all, t);
                 for (size_t extruder_id = 1; extruder_id < num_facets_states; ++extruder_id)
                     if (!sides[bucket][extruder_id].empty())
                         merged[bucket][extruder_id] = tops_all.empty() ? std::move(sides[bucket][extruder_id]) :

src/libslic3r/PrintObject.cpp

@@ -43,6 +43,7 @@
 #include <boost/log/trivial.hpp>
 
 #include <tbb/parallel_for.h>
+#include <tbb/parallel_invoke.h>
 #include <tbb/spin_mutex.h>
 #include <tbb/concurrent_unordered_set.h>
 
@@ -1723,7 +1724,7 @@ void PrintObject::detect_surfaces_type()
                     if (upper_layer) {
                         ExPolygons upper_slices = interface_shells ?
                             diff_ex(layerm_slices_surfaces, upper_layer->m_regions[region_id]->slices.surfaces, ApplySafetyOffset::Yes) :
-                            diff_ex(layerm_slices_surfaces, upper_layer->lslices, ApplySafetyOffset::Yes);
+                            diff_ex_by_piece(layerm_slices_surfaces, to_polygons(upper_layer->lslices), ApplySafetyOffset::Yes);
                         surfaces_append(top, opening_ex(upper_slices, offset), stTop);
                     } else {
                         // if no upper layer, all surfaces of this one are solid
@@ -1749,7 +1750,7 @@ void PrintObject::detect_surfaces_type()
                         surfaces_append(
                             bottom,
                             opening_ex(
-                                diff_ex(layerm_slices_surfaces, lower_layer->lslices, ApplySafetyOffset::Yes),
+                                diff_ex_by_piece(layerm_slices_surfaces, to_polygons(lower_layer->lslices), ApplySafetyOffset::Yes),
                                 offset),
                             surface_type_bottom_other);
                         // if user requested internal shells, we need to identify surfaces
@@ -1780,7 +1781,7 @@ void PrintObject::detect_surfaces_type()
                     // and top surfaces; let's do an intersection to discover them and consider them
                     // as bottom surfaces (to allow for bridge detection)
                     if (! top.empty() && ! bottom.empty()) {
-                        const auto cracks = intersection_ex(top, bottom);
+                        const auto cracks = intersection_ex_by_piece(to_expolygons(top), to_polygons(bottom));
                         if (!cracks.empty()) {
                             if (lower_layer) { // Only detect small cracks for non-first layer, because first layer should always be bottom
                                 const float small_crack_threshold = -layerm->flow(frExternalPerimeter).scaled_width() * 1.5;
@@ -1815,9 +1816,9 @@ void PrintObject::detect_surfaces_type()
                                 }
                             }
 
-                            Polygons top_polygons = to_polygons(std::move(top));
+                            ExPolygons top_expolygons = to_expolygons(std::move(top));
                             top.clear();
-                            surfaces_append(top, diff_ex(top_polygons, bottom), stTop);
+                            surfaces_append(top, diff_ex_by_piece(top_expolygons, to_polygons(bottom)), stTop);
                         }
                     }
 
@@ -1908,7 +1909,7 @@ void PrintObject::detect_surfaces_type()
                     {
                         Polygons topbottom = to_polygons(top);
                         polygons_append(topbottom, to_polygons(bottom));
-                        surfaces_append(surfaces_out, diff_ex(surfaces_prev_expolys, topbottom), stInternal);
+                        surfaces_append(surfaces_out, diff_ex_by_piece(surfaces_prev_expolys, topbottom), stInternal);
                     }
 
                     surfaces_append(surfaces_out, std::move(top));
@@ -2248,10 +2249,10 @@ void PrintObject::discover_vertical_shells()
             // The "ensure vertical wall thickness" feature is not applicable to any of the regions. Quit.
             return;
         BOOST_LOG_TRIVIAL(debug) << "Discovering vertical shells in parallel - start : cache top / bottom";
-        //FIXME Improve the heuristics for a grain size.
+        // One layer per task: on a layer cut through a fine relief the unions below take far longer than elsewhere, and a
-        size_t grain_size = std::max(num_layers / 16, size_t(1));
+        // few such layers next to each other must not end up in one task.
         tbb::parallel_for(
-            tbb::blocked_range<size_t>(0, num_layers, grain_size),
+            tbb::blocked_range<size_t>(0, num_layers, 1),
             [this, &cache_top_botom_regions](const tbb::blocked_range<size_t>& range) {
                 const std::initializer_list<SurfaceType> surfaces_bottom { stBottom, stBottomBridge };
                 const size_t num_regions = this->num_printing_regions();
@@ -2259,56 +2260,66 @@ void PrintObject::discover_vertical_shells()
                     m_print->throw_if_canceled();
                     const Layer                      &layer = *m_layers[idx_layer];
                     DiscoverVerticalShellsCacheEntry &cache = cache_top_botom_regions[idx_layer];
-                    // Simulate single set of perimeters over all merged regions.
+                    const auto top_bottom_expansion = [&layer](size_t region_id) {
-                    float                             perimeter_offset = 0.f;
+                        return float(layer.m_regions[region_id]->flow(frSolidInfill).scaled_spacing()) * top_bottom_expansion_coeff;
-                    float                             perimeter_min_spacing = FLT_MAX;
+                    };
 #ifdef SLIC3R_DEBUG_SLICE_PROCESSING
                     static size_t debug_idx = 0;
                     ++ debug_idx;
 #endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-                    for (size_t region_id = 0; region_id < num_regions; ++ region_id) {
+                    // The top surfaces, the bottom surfaces and the holes are independent of each other.
-                        LayerRegion &layerm               = *layer.m_regions[region_id];
+                    tbb::parallel_invoke(
-                        float        top_bottom_expansion = float(layerm.flow(frSolidInfill).scaled_spacing()) * top_bottom_expansion_coeff;
+                        [&]() {
-                        // Top surfaces.
+                            for (size_t region_id = 0; region_id < num_regions; ++ region_id)
-                        append(cache.top_surfaces, offset(layerm.slices.filter_by_type(stTop), top_bottom_expansion));
+                                append(cache.top_surfaces, offset(layer.m_regions[region_id]->slices.filter_by_type(stTop), top_bottom_expansion(region_id)));
-//                        append(cache.top_surfaces, offset(layerm.fill_surfaces.filter_by_type(stTop), top_bottom_expansion));
+//                              append(cache.top_surfaces, offset(layerm.fill_surfaces.filter_by_type(stTop), top_bottom_expansion));
-                        // Bottom surfaces.
+                            // Save some computing time by reducing the number of polygons.
-                        append(cache.bottom_surfaces, offset(layerm.slices.filter_by_types(surfaces_bottom), top_bottom_expansion));
+                            cache.top_surfaces = union_(cache.top_surfaces);
-//                        append(cache.bottom_surfaces, offset(layerm.fill_surfaces.filter_by_types(surfaces_bottom), top_bottom_expansion));
+                        },
-                        // Calculate the maximum perimeter offset as if the slice was extruded with a single extruder only.
+                        [&]() {
-                        // First find the maxium number of perimeters per region slice.
+                            for (size_t region_id = 0; region_id < num_regions; ++ region_id)
-                        unsigned int perimeters = 0;
+                                append(cache.bottom_surfaces, offset(layer.m_regions[region_id]->slices.filter_by_types(surfaces_bottom), top_bottom_expansion(region_id)));
-                        for (Surface &s : layerm.slices.surfaces)
+//                              append(cache.bottom_surfaces, offset(layerm.fill_surfaces.filter_by_types(surfaces_bottom), top_bottom_expansion));
-                            perimeters = std::max<unsigned int>(perimeters, s.extra_perimeters);
+                            cache.bottom_surfaces = union_(cache.bottom_surfaces);
-                        perimeters += layerm.region().config().wall_loops.value;
+                        },
-                        // Then calculate the infill offset.
+                        [&]() {
-                        if (perimeters > 0) {
+                            // Simulate single set of perimeters over all merged regions.
-                            Flow extflow = layerm.flow(frExternalPerimeter);
+                            float perimeter_offset      = 0.f;
-                            Flow flow    = layerm.flow(frPerimeter);
+                            float perimeter_min_spacing = FLT_MAX;
-                            perimeter_offset = std::max(perimeter_offset,
+                            for (size_t region_id = 0; region_id < num_regions; ++ region_id) {
-                                0.5f * float(extflow.scaled_width() + extflow.scaled_spacing()) + (float(perimeters) - 1.f) * flow.scaled_spacing());
+                                const LayerRegion &layerm = *layer.m_regions[region_id];
-                            perimeter_min_spacing = std::min(perimeter_min_spacing, float(std::min(extflow.scaled_spacing(), flow.scaled_spacing())));
+                                // Calculate the maximum perimeter offset as if the slice was extruded with a single extruder only.
-                        }
+                                // First find the maxium number of perimeters per region slice.
-                        polygons_append(cache.holes, to_polygons(layerm.fill_expolygons));
+                                unsigned int perimeters = 0;
-                    }
+                                for (const Surface &s : layerm.slices.surfaces)
-                    // Save some computing time by reducing the number of polygons.
+                                    perimeters = std::max<unsigned int>(perimeters, s.extra_perimeters);
-                    cache.top_surfaces    = union_(cache.top_surfaces);
+                                perimeters += layerm.region().config().wall_loops.value;
-                    cache.bottom_surfaces = union_(cache.bottom_surfaces);
+                                // Then calculate the infill offset.
-                    // For a multi-material print, simulate perimeter / infill split as if only a single extruder has been used for the whole print.
+                                if (perimeters > 0) {
-                    if (perimeter_offset > 0.) {
+                                    Flow extflow = layerm.flow(frExternalPerimeter);
-                        // The layer.lslices are forced to merge by expanding them first.
+                                    Flow flow    = layerm.flow(frPerimeter);
-                        polygons_append(cache.holes, offset2(layer.lslices, 0.3f * perimeter_min_spacing, - perimeter_offset - 0.3f * perimeter_min_spacing));
+                                    perimeter_offset = std::max(perimeter_offset,
+                                        0.5f * float(extflow.scaled_width() + extflow.scaled_spacing()) + (float(perimeters) - 1.f) * flow.scaled_spacing());
+                                    perimeter_min_spacing = std::min(perimeter_min_spacing, float(std::min(extflow.scaled_spacing(), flow.scaled_spacing())));
+                                }
+                                polygons_append(cache.holes, to_polygons(layerm.fill_expolygons));
+                            }
+                            // For a multi-material print, simulate perimeter / infill split as if only a single extruder has been used for the whole print.
+                            if (perimeter_offset > 0.) {
+                                // The layer.lslices are forced to merge by expanding them first.
+                                polygons_append(cache.holes, offset2(layer.lslices, 0.3f * perimeter_min_spacing, - perimeter_offset - 0.3f * perimeter_min_spacing));
 #ifdef SLIC3R_DEBUG_SLICE_PROCESSING
-                        {
+                                {
-                            Slic3r::SVG svg(debug_out_path("discover_vertical_shells-extra-holes-%d.svg", debug_idx), get_extents(layer.lslices));
+                                    Slic3r::SVG svg(debug_out_path("discover_vertical_shells-extra-holes-%d.svg", debug_idx), get_extents(layer.lslices));
-                            svg.draw(layer.lslices, "blue");
+                                    svg.draw(layer.lslices, "blue");
-                            svg.draw(union_ex(cache.holes), "red");
+                                    svg.draw(union_ex(cache.holes), "red");
-                            svg.draw_outline(union_ex(cache.holes), "black", "blue", scale_(0.05));
+                                    svg.draw_outline(union_ex(cache.holes), "black", "blue", scale_(0.05));
-                            svg.Close();
+                                    svg.Close();
-                        }
+                                }
 #endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
-                    }
+                            }
-                    cache.holes = union_(cache.holes);
+                            cache.holes = union_(cache.holes);
+                        });
                 }
             });
         m_print->throw_if_canceled();
@@ -2606,11 +2617,8 @@ void PrintObject::discover_vertical_shells()
                         Polygons object_volume;
                         Polygons internal_volume;
                         {
-                            Polygons shrinked_bottom_slice = idx_layer > 0 ? to_polygons(m_layers[idx_layer - 1]->lslices) : Polygons{};
+                            if (idx_layer > 0 && idx_layer + 1 < m_layers.size())
-                            Polygons shrinked_upper_slice  = (idx_layer + 1) < m_layers.size() ?
+                                object_volume = to_polygons(intersection_ex_by_piece(m_layers[idx_layer - 1]->lslices, to_polygons(m_layers[idx_layer + 1]->lslices)));
-                                                                 to_polygons(m_layers[idx_layer + 1]->lslices) :
-                                                                 Polygons{};
-                            object_volume = intersection(shrinked_bottom_slice, shrinked_upper_slice);
                             internal_volume = closing(polygonsInternal, SCALED_EPSILON);
                         }
 
@@ -2670,8 +2678,9 @@ void PrintObject::discover_vertical_shells()
 #endif /* SLIC3R_DEBUG_SLICE_PROCESSING */
 
                     // Trim the internal & internalvoid by the shell.
-                    Slic3r::ExPolygons new_internal = diff_ex(layerm->fill_surfaces.filter_by_type(stInternal), regularized_shell);
+                    const Polygons     regularized_shell_polygons = to_polygons(regularized_shell);
-                    Slic3r::ExPolygons new_internal_void = diff_ex(layerm->fill_surfaces.filter_by_type(stInternalVoid), regularized_shell);
+                    Slic3r::ExPolygons new_internal = diff_ex_by_piece(to_expolygons(layerm->fill_surfaces.filter_by_type(stInternal)), regularized_shell_polygons);
+                    Slic3r::ExPolygons new_internal_void = diff_ex_by_piece(to_expolygons(layerm->fill_surfaces.filter_by_type(stInternalVoid)), regularized_shell_polygons);
 
 #ifdef SLIC3R_DEBUG_SLICE_PROCESSING
                     {

tests/libslic3r/test_clipper_utils.cpp

@@ -299,3 +299,47 @@ TEST_CASE("Traversing Clipper PolyTree", "[ClipperUtils]") {
         REQUIRE(count_polys(output) == reference.size());
     }
 }
+
+TEST_CASE("Tiled diff and intersection cover the same area as the plain calls", "[ClipperUtils]") {
+    // A grid of disjoint framed squares, enough of them to be split into several tiles.
+    const int  n    = 40;
+    const coord_t cell = scaled<coord_t>(2.), side = scaled<coord_t>(1.5), frame = scaled<coord_t>(0.3);
+    ExPolygons subject;
+    for (int y = 0; y < n; ++ y)
+        for (int x = 0; x < n; ++ x) {
+            const Point o(x * cell, y * cell);
+            ExPolygon square(Polygon({ o, o + Point(side, 0), o + Point(side, side), o + Point(0, side) }));
+            Polygon hole({ o + Point(frame, frame), o + Point(frame, side - frame), o + Point(side - frame, side - frame), o + Point(side - frame, frame) });
+            square.holes.emplace_back(std::move(hole));
+            subject.emplace_back(std::move(square));
+        }
+    // Clip polygons crossing many squares, one of them large with holes of its own.
+    Polygons clip;
+    const coord_t span = n * cell;
+    for (int i = 0; i < 8; ++ i) {
+        const coord_t y0 = coord_t(i) * span / 8, y1 = y0 + scaled<coord_t>(0.9);
+        clip.emplace_back(Polygon({ Point(- cell, y0), Point(span, y0 + cell * 3), Point(span, y1 + cell * 3), Point(- cell, y1) }));
+    }
+    ExPolygon big(Polygon({ Point(span / 4, span / 4), Point(3 * span / 4, span / 4), Point(3 * span / 4, 3 * span / 4), Point(span / 4, 3 * span / 4) }));
+    for (int i = 0; i < 4; ++ i) {
+        const Point o(span / 4 + scaled<coord_t>(3.1) + i * scaled<coord_t>(9.7), span / 4 + scaled<coord_t>(5.3));
+        big.holes.emplace_back(Polygon({ o, o + Point(0, scaled<coord_t>(20.)), o + Point(scaled<coord_t>(5.), scaled<coord_t>(20.)), o + Point(scaled<coord_t>(5.), 0) }));
+    }
+    polygons_append(clip, to_polygons(big));
+
+    const auto xor_area = [](const ExPolygons &a, const ExPolygons &b) { return area(diff_ex(a, b)) + area(diff_ex(b, a)); };
+    const ApplySafetyOffset safety = GENERATE(ApplySafetyOffset::No, ApplySafetyOffset::Yes);
+    const double tolerance = double(scaled<coord_t>(0.001)) * double(span);
+
+    const ExPolygons diff_plain = diff_ex(subject, clip, safety);
+    const ExPolygons diff_tiled = diff_ex_by_piece(subject, clip, safety);
+    REQUIRE(area(diff_plain) > 0.);
+    CHECK_THAT(area(diff_tiled), Catch::Matchers::WithinRel(area(diff_plain), 1e-9));
+    CHECK(xor_area(diff_tiled, diff_plain) < tolerance);
+
+    const ExPolygons intersection_plain = intersection_ex(subject, clip, safety);
+    const ExPolygons intersection_tiled = intersection_ex_by_piece(subject, clip, safety);
+    REQUIRE(area(intersection_plain) > 0.);
+    CHECK_THAT(area(intersection_tiled), Catch::Matchers::WithinRel(area(intersection_plain), 1e-9));
+    CHECK(xor_area(intersection_tiled, intersection_plain) < tolerance);
+}