NEW: Port of Cura's multi-material interlocking (#5775)

* Init port of Cura's MM interlocking

* Refactor a bit

* Fix crash when bottom surface is multi-color

* Fix crash when boundary avoidance is 0

* Add config

---------

Co-authored-by: zhimin.zeng <zhimin.zeng@bambulab.com>
jira: none
Change-Id: I81cacddf46ad5921a7a2a23fff07cc17addceb6f

src/libslic3r/CMakeLists.txt

@@ -432,6 +432,12 @@ set(lisbslic3r_sources
     FilamentGroup.cpp
     GCode/ToolOrderUtils.hpp
     GCode/ToolOrderUtils.cpp
+    FlushVolPredictor.hpp
+    FlushVolPredictor.cpp
+    Interlocking/InterlockingGenerator.hpp
+    Interlocking/InterlockingGenerator.cpp
+    Interlocking/VoxelUtils.hpp
+    Interlocking/VoxelUtils.cpp
 )
 
 if (APPLE)

src/libslic3r/ClipperUtils.cpp

@@ -830,6 +830,13 @@ Slic3r::ExPolygons union_ex(const Slic3r::ExPolygons& poly1, const Slic3r::ExPol
     return union_ex(expolys);
     }
 
+Slic3r::ExPolygons xor_ex(const Slic3r::ExPolygons &subject, const Slic3r::ExPolygon &clip, ApplySafetyOffset do_safety_offset) {
+    return _clipper_ex(ClipperLib::ctXor, ClipperUtils::ExPolygonsProvider(subject), ClipperUtils::ExPolygonProvider(clip), do_safety_offset);
+}
+Slic3r::ExPolygons xor_ex(const Slic3r::ExPolygons &subject, const Slic3r::ExPolygons &clip, ApplySafetyOffset do_safety_offset) {
+    return _clipper_ex(ClipperLib::ctXor, ClipperUtils::ExPolygonsProvider(subject), ClipperUtils::ExPolygonsProvider(clip), do_safety_offset);
+}
+
 template<typename PathsProvider1, typename PathsProvider2>
 Polylines _clipper_pl_open(ClipperLib::ClipType clipType, PathsProvider1 &&subject, PathsProvider2 &&clip)
 {

src/libslic3r/ClipperUtils.hpp

@@ -563,6 +563,9 @@ Slic3r::ExPolygons union_ex(const Slic3r::ExPolygons& poly1, const Slic3r::ExPol
 ClipperLib::PolyTree union_pt(const Slic3r::Polygons &subject);
 ClipperLib::PolyTree union_pt(const Slic3r::ExPolygons &subject);
 
+Slic3r::ExPolygons xor_ex(const Slic3r::ExPolygons &subject, const Slic3r::ExPolygon &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
+Slic3r::ExPolygons xor_ex(const Slic3r::ExPolygons &subject, const Slic3r::ExPolygons &clip, ApplySafetyOffset do_safety_offset = ApplySafetyOffset::No);
+
 Slic3r::Polygons union_pt_chained_outside_in(const Slic3r::Polygons &subject);
 
 ClipperLib::PolyNodes order_nodes(const ClipperLib::PolyNodes &nodes);

src/libslic3r/Interlocking/InterlockingGenerator.cpp

@@ -0,0 +1,331 @@
+// Copyright (c) 2023 UltiMaker
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "InterlockingGenerator.hpp"
+#include "Layer.hpp"
+
+namespace std {
+template<> struct hash<Slic3r::GridPoint3>
+{
+    size_t operator()(const Slic3r::GridPoint3& pp) const noexcept
+    {
+        static int prime  = 31;
+        int        result = 89;
+        result            = static_cast<int>(result * prime + pp.x());
+        result            = static_cast<int>(result * prime + pp.y());
+        result            = static_cast<int>(result * prime + pp.z());
+        return static_cast<size_t>(result);
+    }
+};
+} // namespace std
+
+
+namespace Slic3r {
+
+void InterlockingGenerator::generate_interlocking_structure(PrintObject* print_object)
+{
+    const auto& config = print_object->config();
+    if (!config.interlocking_beam) {
+        return;
+    }
+
+    const float    rotation           = Geometry::deg2rad(config.interlocking_orientation.value);
+    const coord_t  beam_layer_count   = config.interlocking_beam_layer_count;
+    const int      interface_depth    = config.interlocking_depth;
+    const int      boundary_avoidance = config.interlocking_boundary_avoidance;
+    const coord_t  beam_width         = scaled(config.interlocking_beam_width.value);
+
+    const DilationKernel interface_dilation(GridPoint3(interface_depth, interface_depth, interface_depth), DilationKernel::Type::PRISM);
+
+    const bool           air_filtering = boundary_avoidance > 0;
+    const DilationKernel air_dilation(GridPoint3(boundary_avoidance, boundary_avoidance, boundary_avoidance), DilationKernel::Type::PRISM);
+
+    const coord_t cell_width = beam_width + beam_width;
+    const Vec3crd cell_size(cell_width, cell_width, 2 * beam_layer_count);
+
+    for (size_t region_a_index = 0; region_a_index < print_object->num_printing_regions(); region_a_index++) {
+        const PrintRegion& region_a      = print_object->printing_region(region_a_index);
+        const auto         extruder_nr_a = region_a.extruder(FlowRole::frExternalPerimeter);
+
+        for (size_t region_b_index = region_a_index + 1; region_b_index < print_object->num_printing_regions(); region_b_index++) {
+            const PrintRegion& region_b      = print_object->printing_region(region_b_index);
+            const auto         extruder_nr_b = region_b.extruder(FlowRole::frExternalPerimeter);
+            if (extruder_nr_a == extruder_nr_b) {
+                continue;
+            }
+
+            InterlockingGenerator gen(*print_object, region_a_index, region_b_index, beam_width, boundary_avoidance, rotation, cell_size, beam_layer_count,
+                                      interface_dilation, air_dilation, air_filtering);
+            gen.generateInterlockingStructure();
+        }
+    }
+}
+
+std::pair<ExPolygons, ExPolygons> InterlockingGenerator::growBorderAreasPerpendicular(const ExPolygons& a, const ExPolygons& b, const coord_t& detect) const
+{
+    const coord_t min_line =
+        std::min(print_object.printing_region(region_a_index).flow(print_object, frExternalPerimeter, 0.1).scaled_width(),
+                 print_object.printing_region(region_b_index).flow(print_object, frExternalPerimeter, 0.1).scaled_width());
+
+    const ExPolygons total_shrunk = offset_ex(union_ex(offset_ex(a, min_line), offset_ex(b, min_line)), 2 * -min_line);
+
+    ExPolygons from_border_a = diff_ex(a, total_shrunk);
+    ExPolygons from_border_b = diff_ex(b, total_shrunk);
+
+    ExPolygons temp_a, temp_b;
+    for (coord_t i = 0; i < (detect / min_line) + 2; ++i) {
+        temp_a        = offset_ex(from_border_a, min_line);
+        temp_b        = offset_ex(from_border_b, min_line);
+        from_border_a = diff_ex(temp_a, temp_b);
+        from_border_b = diff_ex(temp_b, temp_a);
+    }
+
+    return {from_border_a, from_border_b};
+}
+
+void InterlockingGenerator::handleThinAreas(const std::unordered_set<GridPoint3>& has_all_meshes) const
+{
+    const coord_t     number_of_beams_detect = boundary_avoidance;
+    const coord_t     number_of_beams_expand = boundary_avoidance - 1;
+    constexpr coord_t rounding_errors        = 5;
+
+    const coord_t max_beam_width = beam_width;
+    const coord_t detect         = (max_beam_width * number_of_beams_detect) + rounding_errors;
+    const coord_t expand         = (max_beam_width * number_of_beams_expand) + rounding_errors;
+    const coord_t close_gaps =
+        std::min(print_object.printing_region(region_a_index).flow(print_object, frExternalPerimeter, 0.1).scaled_width(),
+                 print_object.printing_region(region_b_index).flow(print_object, frExternalPerimeter, 0.1).scaled_width()) / 4;
+
+    // Make an inclusionary polygon, to only actually handle thin areas near actual microstructures (so not in skin for example).
+    std::vector<Polygons> near_interlock_per_layer;
+    near_interlock_per_layer.assign(print_object.layer_count(), Polygons());
+    for (const auto& cell : has_all_meshes) {
+        const auto bottom_corner = vu.toLowerCorner(cell);
+        for (coord_t layer_nr = bottom_corner.z();
+             layer_nr < bottom_corner.z() + cell_size.z() && layer_nr < static_cast<coord_t>(near_interlock_per_layer.size()); ++layer_nr) {
+            near_interlock_per_layer[static_cast<size_t>(layer_nr)].push_back(vu.toPolygon(cell));
+        }
+    }
+    for (auto& near_interlock : near_interlock_per_layer) {
+        near_interlock = offset(union_(closing(near_interlock, rounding_errors)), detect);
+        polygons_rotate(near_interlock, rotation);
+    }
+
+    // Only alter layers when they are present in both meshes, zip should take care if that.
+    for (size_t layer_nr = 0; layer_nr < print_object.layer_count(); layer_nr++){
+        auto       layer   = print_object.get_layer(layer_nr);
+        ExPolygons polys_a = to_expolygons(layer->get_region(region_a_index)->slices.surfaces);
+        ExPolygons polys_b = to_expolygons(layer->get_region(region_b_index)->slices.surfaces);
+
+        const auto [from_border_a, from_border_b] = growBorderAreasPerpendicular(polys_a, polys_b, detect);
+
+        // Get the areas of each mesh that are _not_ thin (large), by performing a morphological open.
+        const ExPolygons large_a = opening_ex(polys_a, detect);
+        const ExPolygons large_b = opening_ex(polys_b, detect);
+
+        // Derive the area that the thin areas need to expand into (so the added areas to the thin strips) from the information we already have.
+        const ExPolygons thin_expansion_a =
+            offset_ex(intersection_ex(intersection_ex(intersection_ex(large_b, offset_ex(diff_ex(polys_a, large_a), expand)),
+                                                      near_interlock_per_layer[layer_nr]),
+                                      from_border_a),
+                      rounding_errors);
+        const ExPolygons thin_expansion_b =
+            offset_ex(intersection_ex(intersection_ex(intersection_ex(large_a, offset_ex(diff_ex(polys_b, large_b), expand)),
+                                                      near_interlock_per_layer[layer_nr]),
+                                      from_border_b),
+                      rounding_errors);
+
+        // Expanded thin areas of the opposing polygon should 'eat into' the larger areas of the polygon,
+        // and conversely, add the expansions to their own thin areas.
+        layer->get_region(region_a_index)->slices.set(closing_ex(diff_ex(union_ex(polys_a, thin_expansion_a), thin_expansion_b), close_gaps), stInternal);
+        layer->get_region(region_b_index)->slices.set(closing_ex(diff_ex(union_ex(polys_b, thin_expansion_b), thin_expansion_a), close_gaps), stInternal);
+    }
+}
+
+void InterlockingGenerator::generateInterlockingStructure() const
+{
+    std::vector<std::unordered_set<GridPoint3>> voxels_per_mesh = getShellVoxels(interface_dilation);
+
+    std::unordered_set<GridPoint3>& has_any_mesh   = voxels_per_mesh[0];
+    std::unordered_set<GridPoint3>& has_all_meshes = voxels_per_mesh[1];
+    has_any_mesh.merge(has_all_meshes); // perform union and intersection simultaneously. Cannibalizes voxels_per_mesh
+
+    if (has_all_meshes.empty()) {
+        return;
+    }
+
+    const std::vector<ExPolygons> layer_regions = computeUnionedVolumeRegions();
+
+    if (air_filtering) {
+        std::unordered_set<GridPoint3> air_cells;
+        addBoundaryCells(layer_regions, air_dilation, air_cells);
+
+        for (const GridPoint3& p : air_cells) {
+            has_all_meshes.erase(p);
+        }
+
+        handleThinAreas(has_all_meshes);
+    }
+
+    applyMicrostructureToOutlines(has_all_meshes, layer_regions);
+}
+
+std::vector<std::unordered_set<GridPoint3>> InterlockingGenerator::getShellVoxels(const DilationKernel& kernel) const
+{
+    std::vector<std::unordered_set<GridPoint3>> voxels_per_mesh(2);
+
+    // mark all cells which contain some boundary
+    for (size_t region_idx = 0; region_idx < 2; region_idx++)
+    {
+        const size_t region = (region_idx == 0) ? region_a_index : region_b_index;
+        std::unordered_set<GridPoint3>& mesh_voxels = voxels_per_mesh[region_idx];
+
+        std::vector<ExPolygons> rotated_polygons_per_layer(print_object.layer_count());
+        for (size_t layer_nr = 0; layer_nr < print_object.layer_count(); layer_nr++)
+        {
+            auto layer = print_object.get_layer(layer_nr);
+            rotated_polygons_per_layer[layer_nr] = to_expolygons(layer->get_region(region)->slices.surfaces);
+            expolygons_rotate(rotated_polygons_per_layer[layer_nr], rotation);
+        }
+
+        addBoundaryCells(rotated_polygons_per_layer, kernel, mesh_voxels);
+    }
+
+    return voxels_per_mesh;
+}
+
+void InterlockingGenerator::addBoundaryCells(const std::vector<ExPolygons>&  layers,
+                                             const DilationKernel&           kernel,
+                                             std::unordered_set<GridPoint3>& cells) const
+{
+    auto voxel_emplacer = [&cells](GridPoint3 p) {
+        if (p.z() < 0) {
+            return true;
+        }
+        cells.emplace(p);
+        return true;
+    };
+
+    for (size_t layer_nr = 0; layer_nr < layers.size(); layer_nr++) {
+        const coord_t z = static_cast<coord_t>(layer_nr);
+        vu.walkDilatedPolygons(layers[layer_nr], z, kernel, voxel_emplacer);
+        ExPolygons skin = layers[layer_nr];
+        if (layer_nr > 0) {
+            skin = xor_ex(skin, layers[layer_nr - 1]);
+        }
+        skin = opening_ex(skin, cell_size.x() / 2.f); // remove superfluous small areas, which would anyway be included because of walkPolygons
+        vu.walkDilatedAreas(skin, z, kernel, voxel_emplacer);
+    }
+}
+
+std::vector<ExPolygons> InterlockingGenerator::computeUnionedVolumeRegions() const
+{
+    const size_t max_layer_count = print_object.layer_count() +
+                                   1; // introduce ghost layer on top for correct skin computation of topmost layer.
+    std::vector<ExPolygons> layer_regions(max_layer_count);
+
+    for (size_t layer_nr = 0; layer_nr < max_layer_count - 1; layer_nr++) {
+        auto& layer_region = layer_regions[static_cast<size_t>(layer_nr)];
+        for (size_t region_idx : {region_a_index, region_b_index}) {
+            auto layer = print_object.get_layer(layer_nr);
+            expolygons_append(layer_region, to_expolygons(layer->get_region(region_idx)->slices.surfaces));
+        }
+        layer_region = closing_ex(layer_region, ignored_gap_); // Morphological close to merge meshes into single volume
+        expolygons_rotate(layer_region, rotation);
+    }
+    return layer_regions;
+}
+
+std::vector<std::vector<ExPolygons>> InterlockingGenerator::generateMicrostructure() const
+{
+    std::vector<std::vector<ExPolygons>> cell_area_per_mesh_per_layer;
+    cell_area_per_mesh_per_layer.resize(2);
+    cell_area_per_mesh_per_layer[0].resize(2);
+    const coord_t beam_w_sum = beam_width + beam_width;
+    const coord_t middle     = cell_size.x() * beam_width / beam_w_sum;
+    const coord_t width[2]   = {middle, cell_size.x() - middle};
+    for (size_t mesh_idx : {0ul, 1ul}) {
+        Point offset(mesh_idx ? middle : 0, 0);
+        Point area_size(width[mesh_idx], cell_size.y());
+
+        Polygon poly;
+        poly.append(offset);
+        poly.append(offset + Point(area_size.x(), 0));
+        poly.append(offset + area_size);
+        poly.append(offset + Point(0, area_size.y()));
+        cell_area_per_mesh_per_layer[0][mesh_idx].emplace_back(poly);
+    }
+    cell_area_per_mesh_per_layer[1] = cell_area_per_mesh_per_layer[0];
+    for (ExPolygons& polys : cell_area_per_mesh_per_layer[1]) {
+        for (ExPolygon& poly : polys) {
+            for (Point& p : poly.contour) {
+                std::swap(p.x(), p.y());
+            }
+        }
+    }
+    return cell_area_per_mesh_per_layer;
+}
+
+void InterlockingGenerator::applyMicrostructureToOutlines(const std::unordered_set<GridPoint3>& cells,
+                                                          const std::vector<ExPolygons>&        layer_regions) const
+{
+    std::vector<std::vector<ExPolygons>> cell_area_per_mesh_per_layer = generateMicrostructure();
+
+    const float  unapply_rotation = -rotation;
+    const size_t max_layer_count  = print_object.layer_count();
+
+    std::vector<ExPolygons> structure_per_layer[2]; // for each mesh the structure on each layer
+
+    // Every `beam_layer_count` number of layers are combined to an interlocking beam layer
+    // to store these we need ceil(max_layer_count / beam_layer_count) of these layers
+    // the formula is rewritten as (max_layer_count + beam_layer_count - 1) / beam_layer_count, so it works for integer division
+    size_t num_interlocking_layers = (max_layer_count + static_cast<size_t>(beam_layer_count) - 1ul) /
+                                     static_cast<size_t>(beam_layer_count);
+    structure_per_layer[0].resize(num_interlocking_layers);
+    structure_per_layer[1].resize(num_interlocking_layers);
+
+    // Only compute cell structure for half the layers, because since our beams are two layers high, every odd layer of the structure will
+    // be the same as the layer below.
+    for (const GridPoint3& grid_loc : cells) {
+        Vec3crd bottom_corner = vu.toLowerCorner(grid_loc);
+        for (size_t mesh_idx = 0; mesh_idx < 2; mesh_idx++) {
+            for (size_t layer_nr = bottom_corner.z(); layer_nr < bottom_corner.z() + cell_size.z() && layer_nr < max_layer_count;
+                 layer_nr += beam_layer_count) {
+                ExPolygons areas_here = cell_area_per_mesh_per_layer[static_cast<size_t>(layer_nr / beam_layer_count) %
+                                                                cell_area_per_mesh_per_layer.size()][mesh_idx];
+                for (auto & here : areas_here) {
+                    here.translate(bottom_corner.x(), bottom_corner.y());
+                }
+                expolygons_append(structure_per_layer[mesh_idx][static_cast<size_t>(layer_nr / beam_layer_count)], areas_here);
+            }
+        }
+    }
+
+    for (size_t mesh_idx = 0; mesh_idx < 2; mesh_idx++) {
+        for (size_t layer_nr = 0; layer_nr < structure_per_layer[mesh_idx].size(); layer_nr++) {
+            ExPolygons& layer_structure = structure_per_layer[mesh_idx][layer_nr];
+            layer_structure = union_ex(layer_structure);
+            expolygons_rotate(layer_structure, unapply_rotation);
+        }
+    }
+
+    for (size_t region_idx = 0; region_idx < 2; region_idx++) {
+        const size_t region = (region_idx == 0) ? region_a_index : region_b_index;
+        for (size_t layer_nr = 0; layer_nr < max_layer_count; layer_nr++) {
+            ExPolygons layer_outlines = layer_regions[layer_nr];
+            expolygons_rotate(layer_outlines, unapply_rotation);
+
+            const ExPolygons areas_here = intersection_ex(structure_per_layer[region_idx][layer_nr / static_cast<size_t>(beam_layer_count)], layer_outlines);
+            const ExPolygons& areas_other = structure_per_layer[!region_idx][layer_nr / static_cast<size_t>(beam_layer_count)];
+
+            auto       layer  = print_object.get_layer(layer_nr);
+            auto&      slices = layer->get_region(region)->slices;
+            ExPolygons polys  = to_expolygons(slices.surfaces);
+            slices.set(union_ex(diff_ex(polys, areas_other), // reduce layer areas inward with beams from other mesh
+                                areas_here)                  // extend layer areas outward with newly added beams
+                       , stInternal);
+        }
+    }
+}
+
+} // namespace Slic3r

src/libslic3r/Interlocking/InterlockingGenerator.hpp

@@ -0,0 +1,172 @@
+// Copyright (c) 2022 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef INTERLOCKING_GENERATOR_HPP
+#define INTERLOCKING_GENERATOR_HPP
+
+#include "../Print.hpp"
+#include "VoxelUtils.hpp"
+
+namespace Slic3r {
+
+/*!
+ * Class for generating an interlocking structure between two adjacent models of a different extruder.
+ *
+ * The structure consists of horizontal beams of the two materials interlaced.
+ * In the z direction the direction of these beams is alternated with 90*.
+ *
+ * Example with two materials # and O
+ * Even beams:      Odd beams:
+ * ######           ##OO##OO
+ * OOOOOO           ##OO##OO
+ * ######           ##OO##OO
+ * OOOOOO           ##OO##OO
+ *
+ * One material of a single cell of the structure looks like this:
+ *                    .-*-.
+ *                .-*       *-.
+ *               |*-.           *-.
+ *               |    *-.           *-.
+ *            .-* *-.     *-.           *-.
+ *        .-*         *-.     *-.       .-*|
+ *    .-*           .-*   *-.     *-.-*    |
+ *   |*-.       .-*     .-*   *-.   |   .-*
+ *   |    *-.-*     .-*           *-|-*
+ *    *-.   |   .-*
+ *        *-|-*
+ *
+ * We set up a voxel grid of (2*beam_w,2*beam_w,2*beam_h) and mark all the voxels which contain both meshes.
+ * We then remove all voxels which also contain air, so that the interlocking pattern will not be visible from the outside.
+ * We then generate and combine the polygons for each voxel and apply those areas to the outlines ofthe meshes.
+ */
+class InterlockingGenerator
+{
+public:
+    /*!
+     * Generate an interlocking structure between each two adjacent meshes.
+     */
+    static void generate_interlocking_structure(PrintObject* print_object);
+
+private:
+    /*!
+     * Generate an interlocking structure between two meshes
+     */
+    void generateInterlockingStructure() const;
+
+    /*!
+     * Private class for storing some variables used in the computation of the interlocking structure between two meshes.
+     * \param region_a_index The first region
+     * \param region_b_index The second region
+     * \param rotation The angle by which to rotate the interlocking pattern
+     * \param cell_size The size of a voxel cell in (coord_t, coord_t, layer_count)
+     * \param beam_layer_count The number of layers for the height of the beams
+     * \param interface_dilation The thicknening kernel for the interface
+     * \param air_dilation The thickening kernel applied to air so that cells near the outside of the model won't be generated
+     * \param air_filtering Whether to fully remove all of the interlocking cells which would be visible on the outside (i.e. touching air).
+     * If no air filtering then those cells will be cut off in the middle of a beam.
+     */
+    InterlockingGenerator(PrintObject&          print_object,
+                          const size_t          region_a_index,
+                          const size_t          region_b_index,
+                          const coord_t         beam_width,
+                          const coord_t         boundary_avoidance,
+                          const float           rotation,
+                          const Vec3crd&        cell_size,
+                          const coord_t         beam_layer_count,
+                          const DilationKernel& interface_dilation,
+                          const DilationKernel& air_dilation,
+                          const bool            air_filtering)
+        : print_object(print_object)
+        , region_a_index(region_a_index)
+        , region_b_index(region_b_index)
+        , beam_width(beam_width)
+        , boundary_avoidance(boundary_avoidance)
+        , vu(cell_size)
+        , rotation(rotation)
+        , cell_size(cell_size)
+        , beam_layer_count(beam_layer_count)
+        , interface_dilation(interface_dilation)
+        , air_dilation(air_dilation)
+        , air_filtering(air_filtering)
+    {}
+
+    /*! Given two polygons, return the parts that border on air, and grow 'perpendicular' up to 'detect' distance.
+     *
+     * \param a The first polygon.
+     * \param b The second polygon.
+     * \param detec The expand distance. (Not equal to offset, but a series of small offsets and differences).
+     * \return A pair of polygons that repressent the 'borders' of a and b, but expanded 'perpendicularly'.
+     */
+    std::pair<ExPolygons, ExPolygons> growBorderAreasPerpendicular(const ExPolygons& a, const ExPolygons& b, const coord_t& detect) const;
+
+    /*! Special handling for thin strips of material.
+     *
+     * Expand the meshes into each other where they need it, namely when a thin strip of material needs to be attached.
+     * \param has_all_meshes Only do this special handling if there's actually microstructure nearby that needs to be adhered to.
+     */
+    void handleThinAreas(const std::unordered_set<GridPoint3>& has_all_meshes) const;
+
+    /*!
+     * Compute the voxels overlapping with the shell of both models.
+     * This includes the walls, but also top/bottom skin.
+     *
+     * \param kernel The dilation kernel to give the returned voxel shell more thickness
+     * \return The shell voxels for mesh a and those for mesh b
+     */
+    std::vector<std::unordered_set<GridPoint3>> getShellVoxels(const DilationKernel& kernel) const;
+
+    /*!
+     * Compute the voxels overlapping with the shell of some layers.
+     * This includes the walls, but also top/bottom skin.
+     *
+     * \param layers The layer outlines for which to compute the shell voxels
+     * \param kernel The dilation kernel to give the returned voxel shell more thickness
+     * \param[out] cells The output cells which elong to the shell
+     */
+    void addBoundaryCells(const std::vector<ExPolygons>& layers, const DilationKernel& kernel, std::unordered_set<GridPoint3>& cells) const;
+
+    /*!
+     * Compute the regions occupied by both models.
+     *
+     * A morphological close is performed so that we don't register small gaps between the two models as being separate.
+     * \return layer_regions The computed layer regions
+     */
+    std::vector<ExPolygons> computeUnionedVolumeRegions() const;
+
+    /*!
+     * Generate the polygons for the beams of a single cell
+     * \return cell_area_per_mesh_per_layer The output polygons for each beam
+     */
+    std::vector<std::vector<ExPolygons>> generateMicrostructure() const;
+
+    /*!
+     * Change the outlines of the meshes with the computed interlocking structure.
+     *
+     * \param cells The cells where we want to apply the interlocking structure.
+     * \param layer_regions The total volume of the two meshes combined (and small gaps closed)
+     */
+    void applyMicrostructureToOutlines(const std::unordered_set<GridPoint3>& cells, const std::vector<ExPolygons>& layer_regions) const;
+
+    static const coord_t ignored_gap_ = 100u; //!< Distance between models to be considered next to each other so that an interlocking structure will be generated there
+
+    PrintObject&  print_object;
+    const size_t  region_a_index;
+    const size_t  region_b_index;
+    const coord_t beam_width;
+    const coord_t boundary_avoidance;
+
+    const VoxelUtils vu;
+
+    const float          rotation;
+    const Vec3crd        cell_size;
+    const coord_t        beam_layer_count;
+    const DilationKernel interface_dilation;
+    const DilationKernel air_dilation;
+    // Whether to fully remove all of the interlocking cells which would be visible on the outside. If no air filtering then those cells
+    // will be cut off midway in a beam.
+    const bool air_filtering;
+};
+
+} // namespace Slic3r
+
+#endif // INTERLOCKING_GENERATOR_HPP

src/libslic3r/Interlocking/VoxelUtils.cpp

@@ -0,0 +1,219 @@
+// Copyright (c) 2022 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#include "VoxelUtils.hpp"
+#include "../Geometry.hpp"
+#include "../Fill/FillRectilinear.hpp"
+#include "../Surface.hpp"
+
+namespace Slic3r
+{
+
+DilationKernel::DilationKernel(GridPoint3 kernel_size, DilationKernel::Type type)
+    : kernel_size_(kernel_size)
+    , type_(type)
+{
+    coord_t mult = kernel_size.x() * kernel_size.y() * kernel_size.z(); // multiplier for division to avoid rounding and to avoid use of floating point numbers
+    relative_cells_.reserve(mult);
+    GridPoint3 half_kernel = kernel_size / 2;
+
+    GridPoint3 start = -half_kernel;
+    GridPoint3 end = kernel_size - half_kernel;
+    for (coord_t x = start.x(); x < end.x(); x++)
+    {
+        for (coord_t y = start.y(); y < end.y(); y++)
+        {
+            for (coord_t z = start.z(); z < end.z(); z++)
+            {
+                GridPoint3 current(x, y, z);
+                if (type != Type::CUBE)
+                {
+                    GridPoint3 limit((x < 0) ? start.x() : end.x() - 1, (y < 0) ? start.y() : end.y() - 1, (z < 0) ? start.z() : end.z() - 1);
+                    if (limit.x() == 0)
+                        limit.x() = 1;
+                    if (limit.y() == 0)
+                        limit.y() = 1;
+                    if (limit.z() == 0)
+                        limit.z() = 1;
+                    const GridPoint3 rel_dists = (mult * current).array() / limit.array();
+                    if ((type == Type::DIAMOND && rel_dists.x() + rel_dists.y() + rel_dists.z() > mult) || (type == Type::PRISM && rel_dists.x() + rel_dists.y() > mult))
+                    {
+                        continue; // don't consider this cell
+                    }
+                }
+                relative_cells_.emplace_back(x, y, z);
+            }
+        }
+    }
+}
+
+bool VoxelUtils::walkLine(Vec3crd start, Vec3crd end, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    Vec3crd diff = end - start;
+
+    const GridPoint3 start_cell = toGridPoint(start);
+    const GridPoint3 end_cell = toGridPoint(end);
+    if (start_cell == end_cell)
+    {
+        return process_cell_func(start_cell);
+    }
+
+    Vec3crd current_cell = start_cell;
+    while (true)
+    {
+        bool continue_ = process_cell_func(current_cell);
+
+        if (! continue_)
+        {
+            return false;
+        }
+
+        int stepping_dim = -1; // dimension in which the line next exits the current cell
+        double percentage_along_line = std::numeric_limits<double>::max();
+        for (int dim = 0; dim < 3; dim++)
+        {
+            if (diff[dim] == 0)
+            {
+                continue;
+            }
+            coord_t crossing_boundary = toLowerCoord(current_cell[dim], dim) + (diff[dim] > 0) * cell_size_[dim];
+            double percentage_along_line_here = (crossing_boundary - start[dim]) / static_cast<double>(diff[dim]);
+            if (percentage_along_line_here < percentage_along_line)
+            {
+                percentage_along_line = percentage_along_line_here;
+                stepping_dim = dim;
+            }
+        }
+        assert(stepping_dim != -1);
+        if (percentage_along_line > 1.0)
+        {
+            // next cell is beyond the end
+            return true;
+        }
+        current_cell[stepping_dim] += (diff[stepping_dim] > 0) ? 1 : -1;
+    }
+    return true;
+}
+
+
+bool VoxelUtils::walkPolygons(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    for (const Polygon& poly : to_polygons(polys))
+    {
+        Point last = poly.back();
+        for (Point p : poly)
+        {
+            bool continue_ = walkLine(Vec3crd(last.x(), last.y(), z), Vec3crd(p.x(), p.y(), z), process_cell_func);
+            if (! continue_)
+            {
+                return false;
+            }
+            last = p;
+        }
+    }
+    return true;
+}
+
+bool VoxelUtils::walkDilatedPolygons(const ExPolygon& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    ExPolygon translated = polys;
+    GridPoint3 k = kernel.kernel_size_;
+    k.x() %= 2;
+    k.y() %= 2;
+    k.z() %= 2;
+    const Vec3crd translation = (Vec3crd(1, 1, 1) - k).array() * cell_size_.array() / 2;
+    if (translation.x() && translation.y())
+    {
+        translated.translate(Point(translation.x(), translation.y()));
+    }
+    return walkPolygons(translated, z + translation.z(), dilate(kernel, process_cell_func));
+}
+
+bool VoxelUtils::walkAreas(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    ExPolygon translated = polys;
+    const Vec3crd translation = -cell_size_ / 2; // offset half a cell so that the dots of spreadDotsArea are centered on the middle of the cell isntead of the lower corners.
+    if (translation.x() && translation.y())
+    {
+        translated.translate(Point(translation.x(), translation.y()));
+    }
+    return _walkAreas(translated, z, process_cell_func);
+}
+
+static Points spreadDotsArea(const ExPolygon& polygons, Point grid_size)
+{
+    std::unique_ptr<Fill> filler(Fill::new_from_type(ipAlignedRectilinear));
+    filler->angle        = Geometry::deg2rad(90.f);
+    filler->spacing      = unscaled(grid_size.x());
+    filler->bounding_box = get_extents(polygons);
+
+    FillParams params;
+    params.density = 1.f;
+    params.anchor_length_max = 0;
+
+    Surface surface(stInternal, polygons);
+    auto    polylines = filler->fill_surface(&surface, params);
+
+    Points result;
+    for (const Polyline& line : polylines) {
+        assert(line.size() == 2);
+        Point a = line[0];
+        Point b = line[1];
+        assert(a.x() == b.x());
+        if (a.y() > b.y()) {
+            std::swap(a, b);
+        }
+        for (coord_t y = a.y() - (a.y() % grid_size.y()) - grid_size.y(); y < b.y(); y += grid_size.y()) {
+            if (y < a.y())
+                continue;
+            result.emplace_back(a.x(), y);
+        }
+    }
+
+    return result;
+}
+
+bool VoxelUtils::_walkAreas(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    Points skin_points = spreadDotsArea(polys, Point(cell_size_.x(), cell_size_.y()));
+    for (Point p : skin_points)
+    {
+        bool continue_ = process_cell_func(toGridPoint(Vec3crd(p.x() + cell_size_.x() / 2, p.y() + cell_size_.y() / 2, z)));
+        if (! continue_)
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+bool VoxelUtils::walkDilatedAreas(const ExPolygon& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    ExPolygon translated = polys;
+    GridPoint3 k = kernel.kernel_size_;
+    k.x() %= 2;
+    k.y() %= 2;
+    k.z() %= 2;
+    const Vec3crd translation = (Vec3crd(1, 1, 1) - k).array() * cell_size_.array() / 2 // offset half a cell when using an even kernel
+                               - cell_size_.array() / 2; // offset half a cell so that the dots of spreadDotsArea are centered on the middle of the cell isntead of the lower corners.
+    if (translation.x() && translation.y())
+    {
+        translated.translate(Point(translation.x(), translation.y()));
+    }
+    return _walkAreas(translated, z + translation.z(), dilate(kernel, process_cell_func));
+}
+
+std::function<bool(GridPoint3)> VoxelUtils::dilate(const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const
+{
+    return [&process_cell_func, &kernel](GridPoint3 loc)
+    {
+        for (const GridPoint3& rel : kernel.relative_cells_)
+        {
+            bool continue_ = process_cell_func(loc + rel);
+            if (! continue_)
+                return false;
+        }
+        return true;
+    };
+}
+} // namespace cura

src/libslic3r/Interlocking/VoxelUtils.hpp

@@ -0,0 +1,212 @@
+// Copyright (c) 2022 Ultimaker B.V.
+// CuraEngine is released under the terms of the AGPLv3 or higher.
+
+#ifndef UTILS_VOXEL_UTILS_H
+#define UTILS_VOXEL_UTILS_H
+
+#include <functional>
+
+#include "../Polygon.hpp"
+#include "../ExPolygon.hpp"
+
+namespace Slic3r
+{
+
+using GridPoint3 = Vec3crd;
+
+/*!
+ * Class for holding the relative positiongs wrt a reference cell on which to perform a dilation.
+ */
+struct DilationKernel
+{
+    /*!
+     * A cubic kernel checks all voxels in a cube around a reference voxel.
+     *  _____
+     * |\ ___\
+     * | |    |
+     *  \|____|
+     *
+     * A diamond kernel uses a manhattan distance to create a diamond shape around a reference voxel.
+     *  /|\
+     * /_|_\
+     * \ | /
+     *  \|/
+     *
+     * A prism kernel is diamond in XY, but extrudes straight in Z around a reference voxel.
+     *   / \
+     *  /   \
+     * |\   /|
+     * | \ / |
+     * |  |  |
+     *  \ | /
+     *   \|/
+     */
+    enum class Type
+    {
+        CUBE,
+        DIAMOND,
+        PRISM
+    };
+    GridPoint3 kernel_size_; //!< Size of the kernel in number of voxel cells
+    Type type_;
+    std::vector<GridPoint3> relative_cells_; //!< All offset positions relative to some reference cell which is to be dilated
+
+    DilationKernel(GridPoint3 kernel_size, Type type);
+};
+
+/*!
+ * Utility class for walking over a 3D voxel grid.
+ *
+ * Contains the math for intersecting voxels with lines, polgons, areas, etc.
+ */
+class VoxelUtils
+{
+public:
+    using grid_coord_t = coord_t;
+
+    Vec3crd cell_size_;
+
+    VoxelUtils(Vec3crd cell_size)
+        : cell_size_(cell_size)
+    {
+    }
+
+    /*!
+     * Process voxels which a line segment crosses.
+     *
+     * \param start Start point of the line
+     * \param end End point of the line
+     * \param process_cell_func Function to perform on each cell the line crosses
+     * \return Whether executing was stopped short as indicated by the \p cell_processing_function
+     */
+    bool walkLine(Vec3crd start, Vec3crd end, const std::function<bool(GridPoint3)>& process_cell_func) const;
+
+    /*!
+     * Process voxels which the line segments of a polygon crosses.
+     *
+     * \warning Voxels may be processed multiple times!
+     *
+     * \param polys The polygons to walk
+     * \param z The height at which the polygons occur
+     * \param process_cell_func Function to perform on each voxel cell
+     * \return Whether executing was stopped short as indicated by the \p cell_processing_function
+     */
+    bool walkPolygons(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const;
+
+    /*!
+     * Process voxels near the line segments of a polygon.
+     * For each voxel the polygon crosses we process each of the offset voxels according to the kernel.
+     *
+     * \warning Voxels may be processed multiple times!
+     *
+     * \param polys The polygons to walk
+     * \param z The height at which the polygons occur
+     * \param process_cell_func Function to perform on each voxel cell
+     * \return Whether executing was stopped short as indicated by the \p cell_processing_function
+     */
+    bool walkDilatedPolygons(const ExPolygon& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const;
+    bool walkDilatedPolygons(const ExPolygons& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const
+    {
+        for (const auto & poly : polys) {
+            if (!walkDilatedPolygons(poly, z, kernel, process_cell_func)) {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+private:
+    /*!
+     * \warning the \p polys is assumed to be translated by half the cell_size in xy already
+     */
+    bool _walkAreas(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const;
+
+public:
+    /*!
+     * Process all voxels inside the area of a polygons object.
+     *
+     * \warning The voxels along the area are not processed. Thin areas might not process any voxels at all.
+     *
+     * \param polys The area to fill
+     * \param z The height at which the polygons occur
+     * \param process_cell_func Function to perform on each voxel cell
+     * \return Whether executing was stopped short as indicated by the \p cell_processing_function
+     */
+    bool walkAreas(const ExPolygon& polys, coord_t z, const std::function<bool(GridPoint3)>& process_cell_func) const;
+
+    /*!
+     * Process all voxels inside the area of a polygons object.
+     * For each voxel inside the polygon we process each of the offset voxels according to the kernel.
+     *
+     * \warning The voxels along the area are not processed. Thin areas might not process any voxels at all.
+     *
+     * \param polys The area to fill
+     * \param z The height at which the polygons occur
+     * \param process_cell_func Function to perform on each voxel cell
+     * \return Whether executing was stopped short as indicated by the \p cell_processing_function
+     */
+    bool walkDilatedAreas(const ExPolygon& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const;
+    bool walkDilatedAreas(const ExPolygons& polys, coord_t z, const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const
+    {
+        for (const auto & poly : polys) {
+            if (!walkDilatedAreas(poly, z, kernel, process_cell_func)) {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+    /*!
+     * Dilate with a kernel.
+     *
+     * Extends the \p process_cell_func, so that for each cell we process nearby cells as well.
+     *
+     * Apply this function to a process_cell_func to create a new process_cell_func which applies the effect to nearby voxels as well.
+     *
+     * \param kernel The offset positions relative to the input of \p process_cell_func
+     * \param process_cell_func Function to perform on each voxel cell
+     */
+    std::function<bool(GridPoint3)> dilate(const DilationKernel& kernel, const std::function<bool(GridPoint3)>& process_cell_func) const;
+
+    GridPoint3 toGridPoint(const Vec3crd& point) const
+    {
+        return GridPoint3(toGridCoord(point.x(), 0), toGridCoord(point.y(), 1), toGridCoord(point.z(), 2));
+    }
+
+    grid_coord_t toGridCoord(const coord_t& coord, const size_t dim) const
+    {
+        assert(dim < 3);
+        return coord / cell_size_[dim] - (coord < 0);
+    }
+
+    Vec3crd toLowerCorner(const GridPoint3& location) const
+    {
+        return Vec3crd(toLowerCoord(location.x(), 0), toLowerCoord(location.y(), 1), toLowerCoord(location.z(), 2));
+    }
+
+    coord_t toLowerCoord(const grid_coord_t& grid_coord, const size_t dim) const
+    {
+        assert(dim < 3);
+        return grid_coord * cell_size_[dim];
+    }
+
+    /*!
+     * Returns a rectangular polygon equal to the cross section of a voxel cell at coordinate \p p
+     */
+    Polygon toPolygon(const GridPoint3 p) const
+    {
+        Polygon ret;
+        Vec3crd c = toLowerCorner(p);
+        ret.append({c.x(), c.y()});
+        ret.append({c.x() + cell_size_.x(), c.y()});
+        ret.append({c.x() + cell_size_.x(), c.y() + cell_size_.y()});
+        ret.append({c.x(), c.y() + cell_size_.y()});
+        return ret;
+    }
+};
+
+} // namespace Slic3r
+
+#endif // UTILS_VOXEL_UTILS_H

src/libslic3r/MultiMaterialSegmentation.cpp

@@ -2309,7 +2309,9 @@ std::vector<std::vector<ExPolygons>> multi_material_segmentation_by_painting(con
     BOOST_LOG_TRIVIAL(debug) << "MM segmentation - layers segmentation in parallel - end";
     throw_on_cancel_callback();
 
-    if (auto max_width = print_object.config().mmu_segmented_region_max_width, interlocking_depth = print_object.config().mmu_segmented_region_interlocking_depth; max_width > 0.f  || interlocking_depth > 0.f) {
+    auto interlocking_beam = print_object.config().interlocking_beam;
+    if (auto max_width = print_object.config().mmu_segmented_region_max_width, interlocking_depth = print_object.config().mmu_segmented_region_interlocking_depth;
+        !interlocking_beam && (max_width > 0.f || interlocking_depth > 0.f)) {
         cut_segmented_layers(input_expolygons, segmented_regions, float(scale_(max_width)), float(scale_(interlocking_depth)), throw_on_cancel_callback);
         throw_on_cancel_callback();
     }

src/libslic3r/Preset.cpp

@@ -900,7 +900,8 @@ static std::vector<std::string> s_Preset_print_options {
     "print_flow_ratio",
     //Orca
     "exclude_object", /*"seam_slope_type",*/ "seam_slope_conditional", "scarf_angle_threshold", /*"seam_slope_start_height", */"seam_slope_entire_loop",/* "seam_slope_min_length",*/
-    "seam_slope_steps", "seam_slope_inner_walls", "role_base_wipe_speed"/*, "seam_slope_gap"*/};
+    "seam_slope_steps", "seam_slope_inner_walls", "role_base_wipe_speed"/*, "seam_slope_gap"*/,
+    "interlocking_beam", "interlocking_orientation", "interlocking_beam_layer_count", "interlocking_depth", "interlocking_boundary_avoidance", "interlocking_beam_width"};
 
 static std::vector<std::string> s_Preset_filament_options {
     /*"filament_colour", */ "default_filament_colour","required_nozzle_HRC","filament_diameter", "filament_type", "filament_soluble", "filament_is_support","filament_scarf_seam_type", "filament_scarf_height", "filament_scarf_gap","filament_scarf_length",

src/libslic3r/PrintConfig.cpp

@@ -2416,6 +2416,56 @@ void PrintConfigDef::init_fff_params()
     def->mode     = comAdvanced;
     def->set_default_value(new ConfigOptionFloat(0.));
 
+    def           = this->add("interlocking_beam", coBool);
+    def->label    = L("Use beam interlocking");
+    def->tooltip  = L("Generate interlocking beam structure at the locations where different filaments touch. This improves the adhesion between filaments, especially models printed in different materials.");
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionBool(false));
+
+    def           = this->add("interlocking_beam_width", coFloat);
+    def->label    = L("Interlocking beam width");
+    def->tooltip  = L("The width of the interlocking structure beams.");
+    def->sidetext = L("mm");
+    def->min      = 0.01;
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionFloat(0.8));
+
+    def           = this->add("interlocking_orientation", coFloat);
+    def->label    = L("Interlocking direction");
+    def->tooltip  = L("Orientation of interlock beams.");
+    def->sidetext = L("°");
+    def->min      = 0;
+    def->max      = 360;
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionFloat(22.5));
+
+    def           = this->add("interlocking_beam_layer_count", coInt);
+    def->label    = L("Interlocking beam layers");
+    def->tooltip  = L("The height of the beams of the interlocking structure, measured in number of layers. Less layers is stronger, but more prone to defects.");
+    def->min      = 1;
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionInt(2));
+
+    def           = this->add("interlocking_depth", coInt);
+    def->label    = L("Interlocking depth");
+    def->tooltip  = L("The distance from the boundary between filaments to generate interlocking structure, measured in cells. Too few cells will result in poor adhesion.");
+    def->min      = 1;
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionInt(2));
+
+    def           = this->add("interlocking_boundary_avoidance", coInt);
+    def->label    = L("Interlocking boundary avoidance");
+    def->tooltip  = L("The distance from the outside of a model where interlocking structures will not be generated, measured in cells.");
+    def->min      = 0;
+    def->category = L("Advanced");
+    def->mode     = comAdvanced;
+    def->set_default_value(new ConfigOptionInt(2));
+
     def = this->add("ironing_type", coEnum);
     def->label = L("Ironing Type");
     def->category = L("Quality");

src/libslic3r/PrintConfig.hpp

@@ -833,6 +833,14 @@ PRINT_CONFIG_CLASS_DEFINE(
     ((ConfigOptionPercent,            wipe_speed))
     ((ConfigOptionBool,               role_base_wipe_speed))
     ((ConfigOptionBool,               precise_z_height)) // BBS
+        
+    ((ConfigOptionBool, interlocking_beam))
+    ((ConfigOptionFloat,interlocking_beam_width))
+    ((ConfigOptionFloat,interlocking_orientation))
+    ((ConfigOptionInt,  interlocking_beam_layer_count))
+    ((ConfigOptionInt,  interlocking_depth))
+    ((ConfigOptionInt,  interlocking_boundary_avoidance))
+
 )
 
 // This object is mapped to Perl as Slic3r::Config::PrintRegion.

src/libslic3r/PrintObject.cpp

@@ -738,7 +738,13 @@ bool PrintObject::invalidate_state_by_config_options(
             || opt_key == "raft_layers"
             || opt_key == "raft_contact_distance"
             || opt_key == "slice_closing_radius"
-            || opt_key == "slicing_mode") {
+            || opt_key == "slicing_mode"
+            || opt_key == "interlocking_beam"
+            || opt_key == "interlocking_orientation"
+            || opt_key == "interlocking_beam_layer_count"
+            || opt_key == "interlocking_depth"
+            || opt_key == "interlocking_boundary_avoidance"
+            || opt_key == "interlocking_beam_width") {
             steps.emplace_back(posSlice);
 		} else if (
                opt_key == "elefant_foot_compensation"

src/libslic3r/PrintObjectSlice.cpp

@@ -4,6 +4,7 @@
 #include "MultiMaterialSegmentation.hpp"
 #include "Print.hpp"
 #include "ClipperUtils.hpp"
+#include "Interlocking/InterlockingGenerator.hpp"
 //BBS
 #include "ShortestPath.hpp"
 
@@ -1070,6 +1071,9 @@ void PrintObject::slice_volumes()
     }
 
 
+    InterlockingGenerator::generate_interlocking_structure(this);
+    m_print->throw_if_canceled();
+
     BOOST_LOG_TRIVIAL(debug) << "Slicing volumes - make_slices in parallel - begin";
     {
         // Compensation value, scaled. Only applying the negative scaling here, as the positive scaling has already been applied during slicing.

src/slic3r/GUI/ConfigManipulation.cpp

@@ -707,6 +707,14 @@ void ConfigManipulation::toggle_print_fff_options(DynamicPrintConfig *config, in
         toggle_field("accel_to_decel_factor", config->opt_bool("accel_to_decel_enable"));
     }
     toggle_line("exclude_object", gcflavor == gcfKlipper);
+
+    bool use_beam_interlocking = config->opt_bool("interlocking_beam");
+    toggle_line("mmu_segmented_region_interlocking_depth", !use_beam_interlocking);
+    toggle_line("interlocking_beam_width", use_beam_interlocking);
+    toggle_line("interlocking_orientation", use_beam_interlocking);
+    toggle_line("interlocking_beam_layer_count", use_beam_interlocking);
+    toggle_line("interlocking_depth", use_beam_interlocking);
+    toggle_line("interlocking_boundary_avoidance", use_beam_interlocking);
 }
 
 void ConfigManipulation::update_print_sla_config(DynamicPrintConfig* config, const bool is_global_config/* = false*/)

src/slic3r/GUI/Tab.cpp

@@ -2213,8 +2213,14 @@ void TabPrint::build()
         optgroup->append_single_option_line("fuzzy_skin_thickness");
 
         optgroup = page->new_optgroup(L("Advanced"), L"advanced");
+        optgroup->append_single_option_line("interlocking_beam");
         // optgroup->append_single_option_line("mmu_segmented_region_max_width");
         optgroup->append_single_option_line("mmu_segmented_region_interlocking_depth");
+        optgroup->append_single_option_line("interlocking_beam_width");
+        optgroup->append_single_option_line("interlocking_orientation");
+        optgroup->append_single_option_line("interlocking_beam_layer_count");
+        optgroup->append_single_option_line("interlocking_depth");
+        optgroup->append_single_option_line("interlocking_boundary_avoidance");
 
         optgroup = page->new_optgroup(L("G-code output"), L"param_gcode");
         optgroup->append_single_option_line("reduce_infill_retraction");