diff --git a/xs/src/libslic3r/ExPolygon.hpp b/xs/src/libslic3r/ExPolygon.hpp
index 095b8f515..4b88e58a1 100644
--- a/xs/src/libslic3r/ExPolygon.hpp
+++ b/xs/src/libslic3r/ExPolygon.hpp
@@ -52,13 +52,14 @@ class ExPolygon
     std::string dump_perl() const;
 };
 
-inline Polygons
-to_polygons(const ExPolygons &expolygons)
+// Count a nuber of polygons stored inside the vector of expolygons.
+// Useful for allocating space for polygons when converting expolygons to polygons.
+inline size_t number_polygons(const ExPolygons &expolys)
 {
-    Polygons pp;
-    for (ExPolygons::const_iterator ex = expolygons.begin(); ex != expolygons.end(); ++ex)
-        append_to(pp, (Polygons)*ex);
-    return pp;
+    size_t n_polygons = 0;
+    for (ExPolygons::const_iterator it = expolys.begin(); it != expolys.end(); ++ it)
+        n_polygons += it->holes.size() + 1;
+    return n_polygons;
 }
 
 inline ExPolygons
@@ -69,8 +70,110 @@ operator+(ExPolygons src1, const ExPolygons &src2) {
 
 std::ostream& operator <<(std::ostream &s, const ExPolygons &expolygons);
 
+
+inline void 
+polygons_append(Polygons &dst, const ExPolygon &src) 
+{ 
+    dst.reserve(dst.size() + src.holes.size() + 1);
+    dst.push_back(src.contour);
+    dst.insert(dst.end(), src.holes.begin(), src.holes.end());
 }
 
+inline void 
+polygons_append(Polygons &dst, const ExPolygons &src) 
+{ 
+    dst.reserve(dst.size() + number_polygons(src));
+    for (ExPolygons::const_iterator it = src.begin(); it != src.end(); ++ it) {
+        dst.push_back(it->contour);
+        dst.insert(dst.end(), it->holes.begin(), it->holes.end());
+    }
+}
+
+inline void 
+polygons_append(Polygons &dst, ExPolygon &&src)
+{ 
+    dst.reserve(dst.size() + src.holes.size() + 1);
+    dst.push_back(std::move(src.contour));
+    std::move(std::begin(src.holes), std::end(src.holes), std::back_inserter(dst));
+    src.holes.clear();
+}
+
+inline void 
+polygons_append(Polygons &dst, ExPolygons &&src)
+{ 
+    dst.reserve(dst.size() + number_polygons(src));
+    for (ExPolygons::iterator it = src.begin(); it != src.end(); ++ it) {
+        dst.push_back(std::move(it->contour));
+        std::move(std::begin(it->holes), std::end(it->holes), std::back_inserter(dst));
+        it->holes.clear();
+    }
+}
+
+inline void 
+expolygons_append(ExPolygons &dst, const ExPolygons &src) 
+{ 
+    dst.insert(dst.end(), src.begin(), src.end());
+}
+
+inline void 
+expolygons_append(ExPolygons &dst, ExPolygons &&src)
+{ 
+    if (dst.empty()) {
+        dst = std::move(src);
+    } else {
+        std::move(std::begin(src), std::end(src), std::back_inserter(dst));
+        src.clear();
+    }
+}
+
+inline Polygons 
+to_polygons(const ExPolygon &src)
+{
+    Polygons polygons;
+    polygons.reserve(src.holes.size() + 1);
+    polygons.push_back(src.contour);
+    polygons.insert(polygons.end(), src.holes.begin(), src.holes.end());
+    return polygons;
+}
+
+inline Polygons 
+to_polygons(const ExPolygons &src)
+{
+    Polygons polygons;
+    polygons.reserve(number_polygons(src));
+    for (ExPolygons::const_iterator it = src.begin(); it != src.end(); ++it) {
+        polygons.push_back(it->contour);
+        polygons.insert(polygons.end(), it->holes.begin(), it->holes.end());
+    }
+    return polygons;
+}
+
+inline Polygons 
+to_polygons(ExPolygon &&src)
+{
+    Polygons polygons;
+    polygons.reserve(src.holes.size() + 1);
+    polygons.push_back(std::move(src.contour));
+    std::move(std::begin(src.holes), std::end(src.holes), std::back_inserter(polygons));
+    src.holes.clear();
+    return polygons;
+}
+
+inline Polygons 
+to_polygons(ExPolygons &&src)
+{
+    Polygons polygons;
+    polygons.reserve(number_polygons(src));
+    for (ExPolygons::iterator it = src.begin(); it != src.end(); ++it) {
+        polygons.push_back(std::move(it->contour));
+        std::move(std::begin(it->holes), std::end(it->holes), std::back_inserter(polygons));
+        it->holes.clear();
+    }
+    return polygons;
+}
+
+} // namespace Slic3r
+
 // start Boost
 #include <boost/polygon/polygon.hpp>
 namespace boost { namespace polygon {
diff --git a/xs/src/libslic3r/Polygon.hpp b/xs/src/libslic3r/Polygon.hpp
index 1534de40c..2e7985d0d 100644
--- a/xs/src/libslic3r/Polygon.hpp
+++ b/xs/src/libslic3r/Polygon.hpp
@@ -21,6 +21,7 @@ class Polygon : public MultiPoint {
     const Point& operator[](Points::size_type idx) const;
     
     Polygon() {};
+
     explicit Polygon(const Points &points): MultiPoint(points) {};
     Point last_point() const;
     virtual Lines lines() const;
@@ -52,6 +53,21 @@ class Polygon : public MultiPoint {
     static Polygon new_scale(const Pointfs& p);
 };
 
+// Append a vector of polygons at the end of another vector of polygons.
+inline void        polygons_append(Polygons &dst, const Polygons &src) { dst.insert(dst.end(), src.begin(), src.end()); }
+
+inline void        polygons_append(Polygons &dst, Polygons &&src) 
+{
+    if (dst.empty()) {
+        dst = std::move(src);
+    } else {
+        std::move(std::begin(src), std::end(src), std::back_inserter(dst));
+        src.clear();
+    }
+}
+
+
+
 inline Polygons
 operator+(Polygons src1, const Polygons &src2) {
     append_to(src1, src2);
diff --git a/xs/src/libslic3r/Print.hpp b/xs/src/libslic3r/Print.hpp
index 612f06030..3e1f868dd 100644
--- a/xs/src/libslic3r/Print.hpp
+++ b/xs/src/libslic3r/Print.hpp
@@ -160,6 +160,12 @@ class PrintObject
     bool has_support_material() const;
     void detect_surfaces_type();
     void process_external_surfaces();
+
+    /// Combine fill surfaces across layers.
+    /// Idempotence of this method is guaranteed by the fact that we don't remove things from
+    /// fill_surfaces but we only turn them into VOID surfaces, thus preserving the boundaries.
+    void combine_infill();
+
     void bridge_over_infill();
     coordf_t adjust_layer_height(coordf_t layer_height) const;
     std::vector<coordf_t> generate_object_layers(coordf_t first_layer_height);
@@ -167,6 +173,18 @@ class PrintObject
     std::vector<ExPolygons> _slice_region(size_t region_id, std::vector<float> z, bool modifier);
     void _make_perimeters();
     void _infill();
+    /// Find all horizontal shells in  this object
+    void discover_horizontal_shells();
+
+    /// Only active if config->infill_only_where_needed. This step trims the sparse infill,
+    /// so it acts as an internal support. It maintains all other infill types intact.
+    /// Here the internal surfaces and perimeters have to be supported by the sparse infill.
+    ///FIXME The surfaces are supported by a sparse infill, but the sparse infill is only as large as the area to support.
+    /// Likely the sparse infill will not be anchored correctly, so it will not work as intended.
+    /// Also one wishes the perimeters to be supported by a full infill.
+    /// Idempotence of this method is guaranteed by the fact that we don't remove things from
+    /// fill_surfaces but we only turn them into VOID surfaces, thus preserving the boundaries.
+    void clip_fill_surfaces();
     
     private:
     Print* _print;
diff --git a/xs/src/libslic3r/PrintObject.cpp b/xs/src/libslic3r/PrintObject.cpp
index 227c32c49..2d3be7c9e 100644
--- a/xs/src/libslic3r/PrintObject.cpp
+++ b/xs/src/libslic3r/PrintObject.cpp
@@ -1095,6 +1095,112 @@ PrintObject::_infill()
     this->state.set_done(posInfill);
 }
 
+void
+PrintObject::combine_infill()
+{
+    // Work on each region separately.
+    for (size_t region_id = 0; region_id < this->print()->regions.size(); ++ region_id) {
+        const PrintRegion *region = this->print()->regions[region_id];
+        const int every = region->config.infill_every_layers.value;
+        if (every < 2 || region->config.fill_density == 0.)
+            continue;
+        // Limit the number of combined layers to the maximum height allowed by this regions' nozzle.
+        //FIXME limit the layer height to max_layer_height
+        double nozzle_diameter = std::min(
+                this->print()->config.nozzle_diameter.get_at(region->config.infill_extruder.value - 1),
+                this->print()->config.nozzle_diameter.get_at(region->config.solid_infill_extruder.value - 1));
+        // define the combinations
+        std::vector<size_t> combine(this->layers.size(), 0);
+        {
+            double current_height = 0.;
+            size_t num_layers = 0;
+            for (size_t layer_idx = 0; layer_idx < this->layers.size(); ++ layer_idx) {
+                const Layer *layer = this->layers[layer_idx];
+                if (layer->id() == 0)
+                    // Skip first print layer (which may not be first layer in array because of raft).
+                    continue;
+                // Check whether the combination of this layer with the lower layers' buffer
+                // would exceed max layer height or max combined layer count.
+                if (current_height + layer->height >= nozzle_diameter + EPSILON || num_layers >= every) {
+                    // Append combination to lower layer.
+                    combine[layer_idx - 1] = num_layers;
+                    current_height = 0.;
+                    num_layers = 0;
+                }
+                current_height += layer->height;
+                ++ num_layers;
+            }
+
+            // Append lower layers (if any) to uppermost layer.
+            combine[this->layers.size() - 1] = num_layers;
+        }
+
+        // loop through layers to which we have assigned layers to combine
+        for (size_t layer_idx = 0; layer_idx < this->layers.size(); ++ layer_idx) {
+            size_t num_layers = combine[layer_idx];
+            if (num_layers <= 1)
+                continue;
+            // Get all the LayerRegion objects to be combined.
+            std::vector<LayerRegion*> layerms;
+            layerms.reserve(num_layers);
+            for (size_t i = layer_idx + 1 - num_layers; i <= layer_idx; ++ i)
+                layerms.emplace_back(this->layers[i]->regions[region_id]);
+            // We need to perform a multi-layer intersection, so let's split it in pairs.
+            // Initialize the intersection with the candidates of the lowest layer.
+            ExPolygons intersection = to_expolygons(layerms.front()->fill_surfaces.filter_by_type(stInternal));
+            // Start looping from the second layer and intersect the current intersection with it.
+            for (size_t i = 1; i < layerms.size(); ++ i)
+                intersection = intersection_ex(
+                        to_polygons(intersection),
+                        to_polygons(layerms[i]->fill_surfaces.filter_by_type(stInternal)),
+                        false);
+            double area_threshold = layerms.front()->infill_area_threshold();
+            if (! intersection.empty() && area_threshold > 0.)
+                intersection.erase(std::remove_if(intersection.begin(), intersection.end(), 
+                            [area_threshold](const ExPolygon &expoly) { return expoly.area() <= area_threshold; }), 
+                        intersection.end());
+            if (intersection.empty())
+                continue;
+            //            Slic3r::debugf "  combining %d %s regions from layers %d-%d\n",
+            //                scalar(@$intersection),
+            //                ($type == S_TYPE_INTERNAL ? 'internal' : 'internal-solid'),
+            //                $layer_idx-($every-1), $layer_idx;
+            // intersection now contains the regions that can be combined across the full amount of layers,
+            // so let's remove those areas from all layers.
+            Polygons intersection_with_clearance;
+            intersection_with_clearance.reserve(intersection.size());
+            float clearance_offset = 
+                0.5f * layerms.back()->flow(frPerimeter).scaled_width() +
+                // Because fill areas for rectilinear and honeycomb are grown 
+                // later to overlap perimeters, we need to counteract that too.
+                ((region->config.fill_pattern == ipRectilinear   ||
+                  region->config.fill_pattern == ipGrid          ||
+                  region->config.fill_pattern == ipHoneycomb) ? 1.5f : 0.5f) * 
+                layerms.back()->flow(frSolidInfill).scaled_width();
+            for (ExPolygon &expoly : intersection)
+                polygons_append(intersection_with_clearance, offset(expoly, clearance_offset));
+            for (LayerRegion *layerm : layerms) {
+                Polygons internal = to_polygons(layerm->fill_surfaces.filter_by_type(stInternal));
+                layerm->fill_surfaces.remove_type(stInternal);
+                layerm->fill_surfaces.append(diff_ex(internal, intersection_with_clearance, false), stInternal);
+                if (layerm == layerms.back()) {
+                    // Apply surfaces back with adjusted depth to the uppermost layer.
+                    Surface templ(stInternal, ExPolygon());
+                    templ.thickness = 0.;
+                    for (LayerRegion *layerm2 : layerms)
+                        templ.thickness += layerm2->layer()->height;
+                    templ.thickness_layers = (unsigned short)layerms.size();
+                    layerm->fill_surfaces.append(intersection, templ);
+                } else {
+                    // Save void surfaces.
+                    layerm->fill_surfaces.append(
+                            intersection_ex(internal, intersection_with_clearance, false),
+                            stInternalVoid);
+                }
+            }
+        }
+    }
+}
 SupportMaterial *
 PrintObject::_support_material()
 {
@@ -1182,4 +1288,77 @@ PrintObject::generate_support_material()
 
 }
 #endif // SLIC3RXS
+void PrintObject::clip_fill_surfaces()
+{
+    if (! this->config.infill_only_where_needed.value ||
+        ! std::any_of(this->print()->regions.begin(), this->print()->regions.end(), 
+            [](const PrintRegion *region) { return region->config.fill_density > 0; }))
+        return;
+
+    // We only want infill under ceilings; this is almost like an
+    // internal support material.
+    // Proceed top-down, skipping the bottom layer.
+    Polygons upper_internal;
+    for (int layer_id = int(this->layers.size()) - 1; layer_id > 0; -- layer_id) {
+        Layer *layer       = this->layers[layer_id];
+        Layer *lower_layer = this->layers[layer_id - 1];
+        // Detect things that we need to support.
+        // Cummulative slices.
+        Polygons slices;
+        for (const ExPolygon &expoly : layer->slices.expolygons)
+            polygons_append(slices, to_polygons(expoly));
+        // Cummulative fill surfaces.
+        Polygons fill_surfaces;
+        // Solid surfaces to be supported.
+        Polygons overhangs;
+        for (const LayerRegion *layerm : layer->regions)
+            for (const Surface &surface : layerm->fill_surfaces.surfaces) {
+                Polygons polygons = to_polygons(surface.expolygon);
+                if (surface.is_solid())
+                    polygons_append(overhangs, polygons);
+                polygons_append(fill_surfaces, std::move(polygons));
+            }
+        Polygons lower_layer_fill_surfaces;
+        Polygons lower_layer_internal_surfaces;
+        for (const LayerRegion *layerm : lower_layer->regions)
+            for (const Surface &surface : layerm->fill_surfaces.surfaces) {
+                Polygons polygons = to_polygons(surface.expolygon);
+                if (surface.surface_type == stInternal || surface.surface_type == stInternalVoid)
+                    polygons_append(lower_layer_internal_surfaces, polygons);
+                polygons_append(lower_layer_fill_surfaces, std::move(polygons));
+            }
+        // We also need to support perimeters when there's at least one full unsupported loop
+        {
+            // Get perimeters area as the difference between slices and fill_surfaces
+            // Only consider the area that is not supported by lower perimeters
+            Polygons perimeters = intersection(diff(slices, fill_surfaces), lower_layer_fill_surfaces);
+            // Only consider perimeter areas that are at least one extrusion width thick.
+            //FIXME Offset2 eats out from both sides, while the perimeters are create outside in.
+            //Should the pw not be half of the current value?
+            float pw = FLT_MAX;
+            for (const LayerRegion *layerm : layer->regions)
+                pw = std::min<float>(pw, layerm->flow(frPerimeter).scaled_width());
+            // Append such thick perimeters to the areas that need support
+            polygons_append(overhangs, offset2(perimeters, -pw, +pw));
+        }
+        // Find new internal infill.
+        polygons_append(overhangs, std::move(upper_internal));
+        upper_internal = intersection(overhangs, lower_layer_internal_surfaces);
+        // Apply new internal infill to regions.
+        for (LayerRegion *layerm : lower_layer->regions) {
+            if (layerm->region()->config.fill_density.value == 0)
+                continue;
+            Polygons internal;
+            for (Surface &surface : layerm->fill_surfaces.surfaces)
+                if (surface.surface_type == stInternal || surface.surface_type == stInternalVoid)
+                    polygons_append(internal, std::move(surface.expolygon));
+            layerm->fill_surfaces.remove_types({ stInternal, stInternalVoid });
+            layerm->fill_surfaces.append(intersection_ex(internal, upper_internal, true), stInternal);
+            layerm->fill_surfaces.append(diff_ex        (internal, upper_internal, true), stInternalVoid);
+            // If there are voids it means that our internal infill is not adjacent to
+            // perimeters. In this case it would be nice to add a loop around infill to
+            // make it more robust and nicer. TODO.
+        }
+    }
+}
 }
diff --git a/xs/src/libslic3r/SurfaceCollection.cpp b/xs/src/libslic3r/SurfaceCollection.cpp
index fb92796ca..b8eb7bb22 100644
--- a/xs/src/libslic3r/SurfaceCollection.cpp
+++ b/xs/src/libslic3r/SurfaceCollection.cpp
@@ -147,5 +147,71 @@ SurfaceCollection::polygons_count() const
         count += 1 + it->expolygon.holes.size();
     return count;
 }
+void
+SurfaceCollection::remove_type(const SurfaceType type)
+{
+    // Use stl remove_if to remove 
+    auto ptr {std::remove_if(surfaces.begin(), surfaces.end(),[type] (Surface& s) { return s.surface_type == type; })};
+    surfaces.erase(ptr, surfaces.end());
+}
+
+void
+SurfaceCollection::remove_types(const SurfaceType *types, size_t ntypes) 
+{
+    for (size_t i = 0; i < ntypes; ++i)
+        this->remove_type(types[i]);
+}
+
+void 
+SurfaceCollection::remove_types(std::initializer_list<SurfaceType> types) {
+    for (const auto& t : types) {
+        this->remove_type(t);
+    }
+}
+
+void
+SurfaceCollection::keep_type(const SurfaceType type)
+{
+    // Use stl remove_if to remove 
+    auto ptr {std::remove_if(surfaces.begin(), surfaces.end(),[type] (Surface& s) { return s.surface_type != type; })};
+    surfaces.erase(ptr, surfaces.end());
+}
+
+void
+SurfaceCollection::keep_types(const SurfaceType *types, size_t ntypes) 
+{
+    for (size_t i = 0; i < ntypes; ++i)
+        this->keep_type(types[i]);
+}
+
+void 
+SurfaceCollection::keep_types(std::initializer_list<SurfaceType> types) {
+    for (const auto& t : types) {
+        this->keep_type(t);
+    }
+}
+/* group surfaces by common properties */
+void
+SurfaceCollection::group(std::vector<SurfacesPtr> *retval)
+{
+    for (Surfaces::iterator it = this->surfaces.begin(); it != this->surfaces.end(); ++it) {
+        // find a group with the same properties
+        SurfacesPtr* group = NULL;
+        for (std::vector<SurfacesPtr>::iterator git = retval->begin(); git != retval->end(); ++git)
+            if (! git->empty() && surfaces_could_merge(*git->front(), *it)) {
+                group = &*git;
+                break;
+            }
+        // if no group with these properties exists, add one
+        if (group == NULL) {
+            retval->resize(retval->size() + 1);
+            group = &retval->back();
+        }
+        // append surface to group
+        group->push_back(&*it);
+    }
+}
+
+
 
 }
diff --git a/xs/src/libslic3r/SurfaceCollection.hpp b/xs/src/libslic3r/SurfaceCollection.hpp
index bb75e0d0d..a12aecec0 100644
--- a/xs/src/libslic3r/SurfaceCollection.hpp
+++ b/xs/src/libslic3r/SurfaceCollection.hpp
@@ -24,6 +24,25 @@ class SurfaceCollection
     SurfacesPtr filter_by_type(SurfaceType type);
     void filter_by_type(SurfaceType type, Polygons* polygons);
 
+    /// deletes all surfaces that match the supplied type.
+    void remove_type(const SurfaceType type);
+
+    void remove_types(std::initializer_list<SurfaceType> types);
+
+    template<int N>
+    void remove_types(std::array<SurfaceType, N> types) {
+        remove_types(types.data(), types.size());
+    }
+    /// group surfaces by common properties
+    void group(std::vector<SurfacesPtr> *retval);
+
+    void keep_type(const SurfaceType type);
+    void keep_types(std::initializer_list<SurfaceType> types);
+    void keep_types(const SurfaceType *types, int ntypes);
+
+    /// deletes all surfaces that match the supplied aggregate of types.
+    void remove_types(const SurfaceType *types, size_t ntypes);
+
     void set(const SurfaceCollection &coll) { surfaces = coll.surfaces; }
     void set(SurfaceCollection &&coll) { surfaces = std::move(coll.surfaces); }
     void set(const ExPolygons &src, SurfaceType surfaceType) { clear(); this->append(src, surfaceType); }
@@ -42,6 +61,10 @@ class SurfaceCollection
     size_t size() const { return this->surfaces.size(); };
     void clear() { this->surfaces.clear(); };
     void erase(size_t i) { this->surfaces.erase(this->surfaces.begin() + i); };
+    Surfaces::iterator begin() { return this->surfaces.begin();}
+    Surfaces::iterator end() { return this->surfaces.end();}
+    Surfaces::const_iterator cbegin() const { return this->surfaces.cbegin();}
+    Surfaces::const_iterator cend() const { return this->surfaces.cend();}
 };
 
 }