Slices algorithm

src/manifold/include/manifold.h

@@ -276,6 +276,7 @@ class Manifold {
    */
   ///@{
   CrossSection Slice(float height = 0) const;
+  std::vector<CrossSection> Slices(float bottomZ, float topZ, int nSlices) const;
   CrossSection Project() const;
   ///@}
 

src/manifold/src/face_op.cpp

@@ -290,6 +290,80 @@ CrossSection Manifold::Impl::Slice(float height) const {
   return CrossSection(polys);
 }
 
+std::vector<CrossSection> Manifold::Impl::Slices(float bottomZ, float topZ, int nSlices) const {
+    std::vector<CrossSection> sections(nSlices);
+    std::vector<float> sliceHeights(nSlices);
+    float step = (topZ - bottomZ) / (nSlices - 1);
+
+    for (int i = 0; i < nSlices; ++i) {
+        sliceHeights[i] = bottomZ + i * step;
+    }
+
+    std::unordered_map<int, std::unordered_set<int>> sliceToTriangles;
+
+    int numTriangles = halfedge_.size() / 3;
+
+    for (int tri = 0; tri < numTriangles; ++tri) {
+        float minZ = std::numeric_limits<float>::infinity();
+        float maxZ = -std::numeric_limits<float>::infinity();
+        for (const int j : {0, 1, 2}) {
+            const float z = vertPos_[halfedge_[3 * tri + j].startVert].z;
+            minZ = std::min(minZ, z);
+            maxZ = std::max(maxZ, z);
+        }
+
+        for (int s = 0; s < nSlices; ++s) {
+            if (minZ <= sliceHeights[s] && maxZ > sliceHeights[s]) {
+                sliceToTriangles[s].insert(tri);
+            }
+        }
+    }
+
+    for (int s = 0; s < nSlices; ++s) {
+        float height = sliceHeights[s];
+        std::unordered_set<int>& tris = sliceToTriangles[s];
+        Polygons polys;
+
+        while (!tris.empty()) {
+            const int startTri = *tris.begin();
+            SimplePolygon poly;
+
+            int k = 0;
+            for (const int j : {0, 1, 2}) {
+                if (vertPos_[halfedge_[3 * startTri + j].startVert].z > height &&
+                    vertPos_[halfedge_[3 * startTri + Next3(j)].startVert].z <= height) {
+                    k = Next3(j);
+                    break;
+                }
+            }
+
+            int tri = startTri;
+            do {
+                tris.erase(tris.find(tri));
+                if (vertPos_[halfedge_[3 * tri + k].endVert].z <= height) {
+                    k = Next3(k);
+                }
+
+                Halfedge up = halfedge_[3 * tri + k];
+                const glm::vec3 below = vertPos_[up.startVert];
+                const glm::vec3 above = vertPos_[up.endVert];
+                const float a = (height - below.z) / (above.z - below.z);
+                poly.push_back(glm::vec2(glm::mix(below, above, a)));
+
+                const int pair = up.pairedHalfedge;
+                tri = pair / 3;
+                k = Next3(pair % 3);
+            } while (tri != startTri);
+
+            polys.push_back(poly);
+        }
+
+        sections[s] = CrossSection(polys);
+    }
+
+    return sections;
+}
+
 CrossSection Manifold::Impl::Project() const {
   const glm::mat3x2 projection = GetAxisAlignedProjection({0, 0, 1});
   auto policy = autoPolicy(halfedge_.size());

src/manifold/src/impl.h

@@ -143,6 +143,7 @@ struct Manifold::Impl {
                             VecView<Halfedge>::IterC end,
                             glm::mat3x2 projection) const;
   CrossSection Slice(float height) const;
+  std::vector<CrossSection> Slices(float bottomZ, float topZ, int nSlices) const;
   CrossSection Project() const;
 
   // edge_op.cu

src/manifold/src/manifold.cpp

@@ -889,6 +889,11 @@ CrossSection Manifold::Slice(float height) const {
   return GetCsgLeafNode().GetImpl()->Slice(height);
 }
 
+std::vector<CrossSection> Manifold::Slices(float bottomZ, float topZ, int nSlices) const {
+  return GetCsgLeafNode().GetImpl()->Slices(bottomZ, topZ, nSlices);
+
+}
+
 /**
  * Returns a cross section representing the projected outline of this object
  * onto the X-Y plane.

test/manifold_test.cpp

@@ -609,6 +609,18 @@ TEST(Manifold, Slice) {
   EXPECT_EQ(top.Area(), 0);
 }
 
+TEST(Manifold, Slices) {
+  Manifold cube = Manifold::Tetrahedron().Scale({10, 10, 10});
+  std::vector<CrossSection> slices = cube.Slices(1.0, 8.0, 8);
+  EXPECT_EQ(slices.size(), 8);
+  EXPECT_GT(slices[0].Area(), 0);
+  for (int i = 0; i < slices.size() - 1; i++) {
+    auto& slice = slices[i];
+    auto& nSlice = slices[i+1];
+    EXPECT_GT(slice.Area(), nSlice.Area());
+  }
+}
+
 TEST(Manifold, MeshRelation) {
   Mesh gyroidMesh = Gyroid();
   MeshGL gyroidMeshGL = WithIndexColors(gyroidMesh);