Add MinkowskiSum and MinkowskiDifference

bindings/python/examples/minkowski.py

@@ -0,0 +1,22 @@
+import numpy as np
+from manifold3d import Manifold
+
+
+def run():
+    small_cube = Manifold.cube([0.1, 0.1, 0.1], True)
+    cube_vertices = small_cube.to_mesh().vert_properties[:, :3]
+    star = Manifold.as_original(small_cube)
+    for offset in [
+        [[0.2, 0.0, 0.0]],
+        [[-0.2, 0.0, 0.0]],
+        [[0.0, 0.2, 0.0]],
+        [[0.0, -0.2, 0.0]],
+        [[0.0, 0.0, 0.2]],
+        [[0.0, 0.0, -0.2]],
+    ]:
+        star += Manifold.hull_points(np.concatenate((cube_vertices, offset), axis=0))
+
+    sphere = Manifold.sphere(0.6, 20)
+    cube = Manifold.cube([1.0, 1.0, 1.0], True)
+    sphereless_cube = cube - sphere
+    return Manifold.minkowski(sphereless_cube, star, False)

bindings/python/manifold3d.cpp

@@ -220,6 +220,9 @@ NB_MODULE(manifold3d, m) {
             return Manifold::Hull(vec);
           },
           "Compute the convex hull enveloping a set of 3d points.")
+      .def_static("minkowski", Manifold::Minkowski, nb::arg("a"), nb::arg("b"),
+                  nb::arg("inset"), nb::arg("useThreading"),
+                  "Compute the minkowski sum of two manifolds.")
       .def(
           "transform",
           [](Manifold &self, nb::ndarray<float, nb::shape<3, 4>> &mat) {

src/manifold/include/manifold.h

@@ -251,6 +251,13 @@ class Manifold {
   static Manifold Hull(const std::vector<glm::vec3>& pts);
   ///@}
 
+  /** @name Minkowski Functions
+   */
+  ///@{
+  static Manifold Minkowski(const Manifold& a, const Manifold& b,
+                            bool inset = false, bool useThreading = false);
+  ///@}
+
   /** @name Testing hooks
    *  These are just for internal testing.
    */

src/manifold/src/face_op.cpp

@@ -316,4 +316,28 @@ CrossSection Manifold::Impl::Project() const {
 
   return CrossSection(polys).Simplify(precision_);
 }
+
+std::vector<int> Manifold::Impl::ReflexFaces(double tolerance) const {
+  std::unordered_set<int> uniqueReflexFaceSet;
+  for (size_t i = 0; i < halfedge_.size(); i++) {
+    Halfedge halfedge = halfedge_[i];
+    int faceA = halfedge.face;
+    int faceB = halfedge_[halfedge.pairedHalfedge].face;
+    glm::dvec3 tangent =
+        glm::cross((glm::dvec3)faceNormal_[faceA],
+                   (glm::dvec3)vertPos_[halfedge_[i].endVert] -
+                       (glm::dvec3)vertPos_[halfedge_[i].startVert]);
+    double tangentProjection =
+        glm::dot((glm::dvec3)faceNormal_[faceB], tangent);
+    //  If we've found a pair of reflex triangles, add them to the set
+    if (tangentProjection > tolerance) {
+      uniqueReflexFaceSet.insert(faceA);
+      uniqueReflexFaceSet.insert(faceB);
+    }
+  }
+  std::vector<int> uniqueFaces;  // Copy to a vector for indexed access
+  uniqueFaces.insert(uniqueFaces.end(), uniqueReflexFaceSet.begin(),
+                     uniqueReflexFaceSet.end());
+  return uniqueFaces;
+}
 }  // namespace manifold

src/manifold/src/impl.h

@@ -121,6 +121,7 @@ struct Manifold::Impl {
                             glm::mat3x2 projection) const;
   CrossSection Slice(float height) const;
   CrossSection Project() const;
+  std::vector<int> ReflexFaces(double tolerance = 1e-8) const;
 
   // edge_op.cu
   void SimplifyTopology();

src/manifold/src/manifold.cpp

@@ -61,6 +61,40 @@ struct UpdateProperties {
   }
 };
 
+struct ComputeTriangleHull {
+  const Manifold* manifold;
+  std::vector<glm::vec3>* vertPos;
+  std::vector<Manifold>* output;
+  void operator()(thrust::tuple<int, glm::ivec3&> inOut) {
+    const int idx = thrust::get<0>(inOut);
+    glm::ivec3& tri = thrust::get<1>(inOut);
+    (*output)[idx] = Manifold::Hull({(*manifold).Translate((*vertPos)[tri.x]),
+                                     (*manifold).Translate((*vertPos)[tri.y]),
+                                     (*manifold).Translate((*vertPos)[tri.z])});
+  }
+};
+
+struct ComputeTriangleTriangleHull {
+  std::vector<glm::vec3>* vertPosAPtr;
+  std::vector<glm::vec3>* vertPosBPtr;
+  std::vector<Manifold>* output;
+  void operator()(thrust::tuple<int, std::pair<glm::ivec3, glm::ivec3>> inOut) {
+    const int idx = thrust::get<0>(inOut);
+    std::pair<glm::ivec3, glm::ivec3> tris = thrust::get<1>(inOut);
+    std::vector<glm::vec3> vertPosA = *vertPosAPtr, vertPosB = *vertPosBPtr;
+    (*output)[idx] =
+        Manifold::Hull({vertPosA[tris.first.x] + vertPosB[tris.second.x],
+                        vertPosA[tris.first.x] + vertPosB[tris.second.y],
+                        vertPosA[tris.first.x] + vertPosB[tris.second.z],
+                        vertPosA[tris.first.y] + vertPosB[tris.second.x],
+                        vertPosA[tris.first.y] + vertPosB[tris.second.y],
+                        vertPosA[tris.first.y] + vertPosB[tris.second.z],
+                        vertPosA[tris.first.z] + vertPosB[tris.second.x],
+                        vertPosA[tris.first.z] + vertPosB[tris.second.y],
+                        vertPosA[tris.first.z] + vertPosB[tris.second.z]});
+  }
+};
+
 Manifold Halfspace(Box bBox, glm::vec3 normal, float originOffset) {
   normal = glm::normalize(normal);
   Manifold cutter =
@@ -859,4 +893,98 @@ Manifold Manifold::Hull() const { return Hull(GetMesh().vertPos); }
 Manifold Manifold::Hull(const std::vector<Manifold>& manifolds) {
   return Compose(manifolds).Hull();
 }
+
+/**
+ * Compute the minkowski sum of two manifolds.
+ *
+ * @param a The first manifold in the sum.
+ * @param b The second manifold in the sum.
+ */
+Manifold Manifold::Minkowski(const Manifold& a, const Manifold& b, bool inset,
+                             bool useThreading) {
+  std::vector<Manifold> composedHulls({a});
+  bool aConvex = a.GetCsgLeafNode().GetImpl()->ReflexFaces().size() == 0;
+  bool bConvex = b.GetCsgLeafNode().GetImpl()->ReflexFaces().size() == 0;
+
+  // If the convex manifold was supplied first, swap them!
+  Manifold aM = a, bM = b;
+  if (aConvex && !bConvex) {
+    aM = b;
+    bM = a;
+    aConvex = !aConvex;
+    bConvex = !bConvex;
+  }
+
+  manifold::Mesh aMesh = aM.GetMesh();
+
+  // Convex-Convex Minkowski: Very Fast
+  if (aConvex && bConvex) {
+    std::vector<Manifold> simpleHull;
+    for (glm::vec3 vertex : aMesh.vertPos) {
+      simpleHull.push_back(bM.Translate(vertex));
+    }
+    composedHulls.push_back(Manifold::Hull(simpleHull));
+    // Convex - Non-Convex Minkowski: Slower
+  } else if (!aConvex && bConvex) {
+    if (useThreading) {
+      composedHulls.resize(aMesh.triVerts.size() + 1);
+      thrust::for_each_n(
+          thrust::host, zip(countAt(1), aMesh.triVerts.begin()),
+          aMesh.triVerts.size(),
+          ComputeTriangleHull({&bM, &aMesh.vertPos, &composedHulls}));
+    } else {
+      std::vector<std::vector<Manifold>> composedParts;
+      for (glm::ivec3 vertexIndices : aMesh.triVerts) {
+        composedParts.push_back({bM.Translate(aMesh.vertPos[vertexIndices.x]),
+                                 bM.Translate(aMesh.vertPos[vertexIndices.y]),
+                                 bM.Translate(aMesh.vertPos[vertexIndices.z])});
+      }
+      std::vector<Manifold> newHulls;
+      newHulls.reserve(composedParts.size());
+      newHulls.resize(composedParts.size());
+      thrust::for_each_n(
+          thrust::host, zip(composedParts.begin(), newHulls.begin()),
+          composedParts.size(),
+          [](thrust::tuple<std::vector<Manifold>, Manifold&> inOut) {
+            thrust::get<1>(inOut) = Manifold::Hull(thrust::get<0>(inOut));
+          });
+      composedHulls.insert(composedHulls.end(), newHulls.begin(),
+                           newHulls.end());
+    }
+    // Non-Convex - Non-Convex Minkowski: Very Slow
+  } else if (!aConvex && !bConvex) {
+    manifold::Mesh bMesh = bM.GetMesh();
+    if (useThreading) {
+      std::vector<std::pair<glm::ivec3, glm::ivec3>> trianglePairs;
+      for (glm::ivec3 aVertexIndices : aMesh.triVerts) {
+        for (glm::ivec3 bVertexIndices : bMesh.triVerts) {
+          trianglePairs.push_back({aVertexIndices, bVertexIndices});
+        }
+      }
+      composedHulls.resize(trianglePairs.size() + 1);
+      thrust::for_each_n(thrust::host, zip(countAt(1), trianglePairs.begin()),
+                         trianglePairs.size(),
+                         ComputeTriangleTriangleHull(
+                             {&aMesh.vertPos, &bMesh.vertPos, &composedHulls}));
+    } else {
+      for (glm::ivec3 aIndices : aMesh.triVerts) {
+        for (glm::ivec3 bIndices : bMesh.triVerts) {
+          composedHulls.push_back(Manifold::Hull(
+              {aMesh.vertPos[aIndices.x] + bMesh.vertPos[bIndices.x],
+               aMesh.vertPos[aIndices.x] + bMesh.vertPos[bIndices.y],
+               aMesh.vertPos[aIndices.x] + bMesh.vertPos[bIndices.z],
+               aMesh.vertPos[aIndices.y] + bMesh.vertPos[bIndices.x],
+               aMesh.vertPos[aIndices.y] + bMesh.vertPos[bIndices.y],
+               aMesh.vertPos[aIndices.y] + bMesh.vertPos[bIndices.z],
+               aMesh.vertPos[aIndices.z] + bMesh.vertPos[bIndices.x],
+               aMesh.vertPos[aIndices.z] + bMesh.vertPos[bIndices.y],
+               aMesh.vertPos[aIndices.z] + bMesh.vertPos[bIndices.z]}));
+        }
+      }
+    }
+  }
+  return Manifold::BatchBoolean(composedHulls, inset
+                                                   ? manifold::OpType::Subtract
+                                                   : manifold::OpType::Add);
+}
 }  // namespace manifold

test/manifold_test.cpp

@@ -784,3 +784,74 @@ TEST(Manifold, EmptyHull) {
       {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}};
   EXPECT_TRUE(Manifold::Hull(coplanar).IsEmpty());
 }
+
+TEST(Manifold, ConvexConvexMinkowski) {
+  Manifold sphere = Manifold::Sphere(0.1, 20);
+  Manifold cube = Manifold::Cube();
+  Manifold sum = Manifold::Minkowski(cube, sphere);
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.6966598f);
+}
+
+TEST(Manifold, ConvexConvexMinkowskiThreaded) {
+  Manifold sphere = Manifold::Sphere(0.1, 20);
+  Manifold cube = Manifold::Cube();
+  Manifold sum = Manifold::Minkowski(cube, sphere, true);
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.6966598f);
+}
+
+TEST(Manifold, NonConvexConvexMinkowski) {
+  Manifold sphere = Manifold::Sphere(0.6, 20);
+  Manifold cube = Manifold::Cube({1.0, 1.0, 1.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum = Manifold::Minkowski(nonConvex, Manifold::Sphere(0.1, 20));
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.0686193f);
+}
+
+TEST(Manifold, NonConvexConvexMinkowskiThreaded) {
+  Manifold sphere = Manifold::Sphere(0.6, 20);
+  Manifold cube = Manifold::Cube({1.0, 1.0, 1.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum =
+      Manifold::Minkowski(nonConvex, Manifold::Sphere(0.1, 20), true);
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.0686193f);
+}
+
+TEST(Manifold, NonConvexNonConvexMinkowski) {
+  ManifoldParams().deterministic = true;
+  Manifold star = Manifold::Cube({0.1, 0.1, 0.1}, true);
+  std::vector<glm::vec3> verts = star.GetMesh().vertPos;
+  for (glm::vec3 point :
+       {glm::vec3(0.2, 0.0, 0.0), glm::vec3(-0.2, 0.0, 0.0),
+        glm::vec3(0.0, 0.2, 0.0), glm::vec3(0.0, -0.2, 0.0),
+        glm::vec3(0.0, 0.0, 0.2), glm::vec3(0.0, 0.0, -0.2)}) {
+    verts.push_back(point);
+    star += Manifold::Hull(verts);
+    verts.pop_back();
+  }
+
+  Manifold sphere = Manifold::Sphere(0.6, 20);
+  Manifold cube = Manifold::Cube({1.0, 1.0, 1.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum = Manifold::Minkowski(nonConvex, star);
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.643248);
+}
+
+TEST(Manifold, NonConvexNonConvexMinkowskiThreaded) {
+  ManifoldParams().deterministic = true;
+  Manifold star = Manifold::Cube({0.1, 0.1, 0.1}, true);
+  std::vector<glm::vec3> verts = star.GetMesh().vertPos;
+  for (glm::vec3 point :
+       {glm::vec3(0.2, 0.0, 0.0), glm::vec3(-0.2, 0.0, 0.0),
+        glm::vec3(0.0, 0.2, 0.0), glm::vec3(0.0, -0.2, 0.0),
+        glm::vec3(0.0, 0.0, 0.2), glm::vec3(0.0, 0.0, -0.2)}) {
+    verts.push_back(point);
+    star += Manifold::Hull(verts);
+    verts.pop_back();
+  }
+
+  Manifold sphere = Manifold::Sphere(0.6, 20);
+  Manifold cube = Manifold::Cube({1.0, 1.0, 1.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum = Manifold::Minkowski(nonConvex, star, true);
+  EXPECT_FLOAT_EQ(sum.GetProperties().volume, 1.643248);
+}