diff --git a/xs/Build.PL b/xs/Build.PL
index 5c58747dd..2a7dc4fa1 100644
--- a/xs/Build.PL
+++ b/xs/Build.PL
@@ -24,7 +24,7 @@ my $build = Module::Build::WithXSpp->new(
     # _GLIBCXX_USE_C99 : to get the long long type for g++
     # HAS_BOOL         : stops Perl/lib/CORE/handy.h from doing "#  define bool char" for MSVC
     # NOGDI            : prevents inclusion of wingdi.h which defines functions Polygon() and Polyline() in global namespace
-    extra_compiler_flags => [qw(-D_GLIBCXX_USE_C99 -DHAS_BOOL -DNOGDI)],
+    extra_compiler_flags => [qw(-D_GLIBCXX_USE_C99 -DHAS_BOOL -DNOGDI), ($ENV{SLIC3R_DEBUG} ? ' -DSLIC3R_DEBUG' : '')],
     
     # Provides extra C typemaps that are auto-merged
     extra_typemap_modules => {
diff --git a/xs/src/TriangleMesh.cpp b/xs/src/TriangleMesh.cpp
index 01df14774..165c32af4 100644
--- a/xs/src/TriangleMesh.cpp
+++ b/xs/src/TriangleMesh.cpp
@@ -1,4 +1,6 @@
 #include "TriangleMesh.hpp"
+#include <algorithm>
+#include <math.h>
 
 namespace Slic3r {
 
@@ -130,4 +132,94 @@ void TriangleMesh::rotate(double angle, Point* center)
     this->translate(+center->x, +center->y, 0);
 }
 
+std::vector<Polygons>*
+TriangleMesh::slice(const std::vector<double> &z)
+{
+    /*
+       This method gets called with a list of Z coordinates and outputs
+       a vector pointer having the same number of items as the original list.
+       Each item is a vector of polygons created by slicing our mesh at the 
+       given heights.
+       
+       This method should basically combine the behavior of the existing
+       Perl methods defined in lib/Slic3r/TriangleMesh.pm:
+       
+       - analyze(): this creates the 'facets_edges' and the 'edges_facets'
+            tables (we don't need the 'edges' table)
+       
+       - slice_facet(): this has to be done for each facet. It generates 
+            intersection lines with each plane identified by the Z list.
+            The get_layer_range() binary search used to identify the Z range
+            of the facet is already ported to C++ (see Object.xsp)
+       
+       - make_loops(): this has to be done for each layer. It creates polygons
+            from the lines generated by the previous step.
+        
+        At the end, we free the tables generated by analyze() as we don't 
+        need them anymore.
+        FUTURE: parallelize slice_facet() and make_loops()
+    */
+    
+    for (int facet_idx = 0; facet_idx < this->stl.stats.number_of_facets; facet_idx++) {
+        stl_facet facet = this->stl.facet_start[facet_idx];  // this is a copy
+        
+        /* reorder vertices so that the first one is the one with lowest Z
+           this is needed to get all intersection lines in a consistent order
+           (external on the right of the line) */
+        float min_z;
+        if (facet.vertex[1].z < facet.vertex[0].z && facet.vertex[1].z < facet.vertex[2].z) {
+            // vertex 1 has lowest Z
+            min_z = facet.vertex[1].z;
+            stl_vertex v0 = facet.vertex[0];
+            facet.vertex[0] = facet.vertex[1];
+            facet.vertex[1] = facet.vertex[2];
+            facet.vertex[2] = v0;
+        } else if (facet.vertex[2].z < facet.vertex[0].z && facet.vertex[2].z < facet.vertex[1].z) {
+            // vertex 2 has lowest Z
+            min_z = facet.vertex[2].z;
+            stl_vertex v0 = facet.vertex[0];
+            facet.vertex[0] = facet.vertex[2];
+            facet.vertex[2] = facet.vertex[1];
+            facet.vertex[1] = v0;
+        } else {
+            min_z = facet.vertex[0].z;
+        }
+        float max_z = fmaxf(facet.vertex[0].z, fmaxf(facet.vertex[1].z, facet.vertex[2].z));
+        
+        #ifdef SLIC3R_DEBUG
+        printf("\n==> FACET %d (%f,%f,%f - %f,%f,%f - %f,%f,%f):\n", facet_idx,
+            facet.vertex[0].x, facet.vertex[0].y, facet.vertex[0].z,
+            facet.vertex[1].x, facet.vertex[1].y, facet.vertex[1].z,
+            facet.vertex[2].x, facet.vertex[2].y, facet.vertex[2].z);
+        printf("z: min = %.2f, max = %.2f\n", min_z, max_z);
+        #endif
+        
+        if (min_z == max_z) {
+            #ifdef SLIC3R_DEBUG
+            printf("Facet is horizontal; ignoring\n");
+            #endif
+            continue;
+        }
+        
+        std::vector<double>::const_iterator min_layer, max_layer;
+        min_layer = std::lower_bound(z.begin(), z.end(), min_z); // first layer whose slice_z is >= min_z
+        max_layer = std::upper_bound(z.begin() + (min_layer - z.begin()), z.end(), max_z) - 1; // last layer whose slice_z is <= max_z
+        #ifdef SLIC3R_DEBUG
+        printf("layers: min = %d, max = %d\n", (int)(min_layer - z.begin()), (int)(max_layer - z.begin()));
+        #endif
+        
+        for (std::vector<double>::const_iterator it = z.begin(); it != z.end(); ++it) {
+            
+        }
+    }
+    
+    std::vector<Polygons>* layers = new std::vector<Polygons>(z.size());
+    
+    // ...
+    // add a Polygon p to layer n:
+    // (*layers)[n].push_back(p);
+    
+    return layers;
+}
+
 }
diff --git a/xs/src/TriangleMesh.hpp b/xs/src/TriangleMesh.hpp
index da3d10451..e6ed3fa0c 100644
--- a/xs/src/TriangleMesh.hpp
+++ b/xs/src/TriangleMesh.hpp
@@ -3,7 +3,9 @@
 
 #include <myinit.h>
 #include <admesh/stl.h>
+#include <vector>
 #include "Point.hpp"
+#include "Polygon.hpp"
 
 namespace Slic3r {
 
@@ -20,6 +22,7 @@ class TriangleMesh
     void translate(float x, float y, float z);
     void align_to_origin();
     void rotate(double angle, Point* center);
+    std::vector<Polygons>* slice(const std::vector<double> &z);
     stl_file stl;
 };
 
diff --git a/xs/t/01_trianglemesh.t b/xs/t/01_trianglemesh.t
index 2cdfe01ef..1948d0be5 100644
--- a/xs/t/01_trianglemesh.t
+++ b/xs/t/01_trianglemesh.t
@@ -44,4 +44,11 @@ my $cube = {
     ok abs($m->size->[0] - sqrt(2)*40) < 1E-4, 'rotate';
 }
 
+{
+    my $m = Slic3r::TriangleMesh::XS->new;
+    $m->ReadFromPerl($cube->{vertices}, $cube->{facets});
+    $m->Repair;
+    my $result = $m->slice([2,4,8,6,8,10,12,14,16,18]);
+}
+
 __END__
diff --git a/xs/xsp/TriangleMesh.xsp b/xs/xsp/TriangleMesh.xsp
index 2c88be3d9..97f66bfae 100644
--- a/xs/xsp/TriangleMesh.xsp
+++ b/xs/xsp/TriangleMesh.xsp
@@ -92,6 +92,28 @@ TriangleMesh::size()
     OUTPUT:
         RETVAL
 
+SV*
+TriangleMesh::slice(z)
+    std::vector<double>* z
+    CODE:
+        std::vector<Polygons>* layers = THIS->slice(*z);
+        
+        AV* layers_av = newAV();
+        av_extend(layers_av, layers->size()-1);
+        for (unsigned int i = 0; i < layers->size(); i++) {
+            AV* polygons_av = newAV();
+            av_extend(polygons_av, (*layers)[i].size()-1);
+            unsigned int j = 0;
+            for (Polygons::iterator it = (*layers)[i].begin(); it != (*layers)[i].end(); ++it) {
+                av_store(polygons_av, j++, (*it).to_SV_clone_ref());
+            }
+            av_store(layers_av, i, newRV_noinc((SV*)polygons_av));
+        }
+        delete layers;
+        RETVAL = (SV*)newRV_noinc((SV*)layers_av);
+    OUTPUT:
+        RETVAL
+
 %}
 };
 
diff --git a/xs/xsp/typemap.xspt b/xs/xsp/typemap.xspt
index eeb4f093b..ee1da30e0 100644
--- a/xs/xsp/typemap.xspt
+++ b/xs/xsp/typemap.xspt
@@ -1,6 +1,7 @@
 %typemap{bool}{simple};
 %typemap{std::string};
 %typemap{std::vector<unsigned int>*};
+%typemap{std::vector<double>*};
 %typemap{SV*};
 %typemap{AV*};
 %typemap{Point*};