Add MinkowskiSum and MinkowskiDifference

bindings/c/include/manifoldc.h

@@ -110,6 +110,10 @@ ManifoldManifoldPair manifold_split_by_plane(void *mem_first, void *mem_second,
 ManifoldManifold *manifold_trim_by_plane(void *mem, ManifoldManifold *m,
                                          float normal_x, float normal_y,
                                          float normal_z, float offset);
+ManifoldManifold *manifold_minkowski_sum(void *mem, ManifoldManifold *a,
+                                         ManifoldManifold *b);
+ManifoldManifold *manifold_minkowski_difference(void *mem, ManifoldManifold *a,
+                                                ManifoldManifold *b);
 
 // 3D to 2D
 

bindings/c/manifoldc.cpp

@@ -202,6 +202,18 @@ ManifoldManifold *manifold_trim_by_plane(void *mem, ManifoldManifold *m,
   return to_c(new (mem) Manifold(trimmed));
 }
 
+ManifoldManifold *manifold_minkowski_sum(void *mem, ManifoldManifold *a,
+                                         ManifoldManifold *b) {
+  auto m = (*from_c(a)).MinkowskiSum(*from_c(b));
+  return to_c(new (mem) Manifold(m));
+}
+
+ManifoldManifold *manifold_minkowski_difference(void *mem, ManifoldManifold *a,
+                                                ManifoldManifold *b) {
+  auto m = (*from_c(a)).MinkowskiDifference(*from_c(b));
+  return to_c(new (mem) Manifold(m));
+}
+
 ManifoldCrossSection *manifold_slice(void *mem, ManifoldManifold *m,
                                      float height) {
   auto poly = from_c(m)->Slice(height);

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 sphereless_cube.minkowski_sum(star)

bindings/python/manifold3d.cpp

@@ -518,6 +518,19 @@ NB_MODULE(manifold3d, m) {
           "vector from the plane.\n"
           ":param origin_offset: The distance of the plane from the origin in "
           "the direction of the normal vector.")
+      .def("minkowski_sum", &Manifold::MinkowskiSum, nb::arg("other"),
+           "Compute the minkowski sum of this manifold with another."
+           "This corresponds to the morphological dilation of the manifold."
+           "\n\n"
+           ":param other: The other manifold to minkowski sum to this one.")
+      .def("minkowski_difference", &Manifold::MinkowskiDifference,
+           nb::arg("other"),
+           "Subtract the sweep of the other manifold across this manifold's "
+           "surface."
+           "This corresponds to the morphological erosion of the manifold."
+           "\n\n"
+           ":param other: The other manifold to minkowski subtract from this "
+           "one.")
       .def("slice", &Manifold::Slice, nb::arg("height"),
            "Returns the cross section of this object parallel to the X-Y plane "
            "at the specified height. Using a height equal to the bottom of the "

bindings/wasm/bindings.cpp

@@ -138,6 +138,8 @@ EMSCRIPTEN_BINDINGS(whatever) {
       .function("_Split", &man_js::Split)
       .function("_SplitByPlane", &man_js::SplitByPlane)
       .function("_TrimByPlane", &Manifold::TrimByPlane)
+      .function("minkowskiSum", &Manifold::MinkowskiSum)
+      .function("minkowskiDifference", &Manifold::MinkowskiDifference)
       .function("slice", &Manifold::Slice)
       .function("project", &Manifold::Project)
       .function("hull", select_overload<Manifold() const>(&Manifold::Hull))

bindings/wasm/examples/worker.ts

@@ -59,9 +59,13 @@ const manifoldStaticFunctions = [
 ];
 // manifold member functions (that return a new manifold)
 const manifoldMemberFunctions = [
-  'add', 'subtract', 'intersect', 'decompose', 'warp', 'transform', 'translate',
-  'rotate', 'scale', 'mirror', 'refine', 'setProperties', 'asOriginal',
-  'trimByPlane', 'split', 'splitByPlane', 'slice', 'project', 'hull'
+  'add',          'subtract',     'intersect',
+  'decompose',    'warp',         'transform',
+  'translate',    'rotate',       'scale',
+  'mirror',       'refine',       'setProperties',
+  'asOriginal',   'trimByPlane',  'split',
+  'splitByPlane', 'slice',        'project',
+  'hull',         'minkowskiSum', 'minkowskiDifference'
 ];
 // CrossSection static methods (that return a new cross-section)
 const crossSectionStaticFunctions = [

bindings/wasm/manifold-encapsulated-types.d.ts

@@ -685,6 +685,22 @@ export class Manifold {
    */
   trimByPlane(normal: Vec3, originOffset: number): Manifold;
 
+  /**
+   * Compute the minkowski sum of this manifold with another.
+   * This corresponds to the morphological dilation of the manifold.
+   *
+   * @param other The other manifold to minkowski sum to this one.
+   */
+  minkowskiSum(other: Manifold): Manifold;
+
+  /**
+   * Subtract the sweep of the other manifold across this manifold's surface.
+   * This corresponds to the morphological erosion of the manifold.
+   *
+   * @param other The other manifold to minkowski subtract from this one.
+   */
+  minkowskiDifference(other: Manifold): Manifold;
+
   /**
    * Returns the cross section of this object parallel to the X-Y plane at the
    * specified height. Using a height equal to the bottom

src/manifold/include/manifold.h

@@ -234,6 +234,8 @@ class Manifold {
   std::pair<Manifold, Manifold> SplitByPlane(glm::vec3 normal,
                                              float originOffset) const;
   Manifold TrimByPlane(glm::vec3 normal, float originOffset) const;
+  Manifold MinkowskiSum(const Manifold&) const;
+  Manifold MinkowskiDifference(const Manifold&) const;
   ///@}
 
   /** @name 2D from 3D
@@ -269,6 +271,8 @@ class Manifold {
   mutable std::shared_ptr<CsgNode> pNode_;
 
   CsgLeafNode& GetCsgLeafNode() const;
+
+  Manifold Minkowski(const Manifold&, bool inset = false) const;
 };
 /** @} */
 }  // namespace manifold

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

@@ -793,6 +793,96 @@
   return *this ^ Halfspace(BoundingBox(), normal, originOffset);
 }
 
+/**
+ * Compute the minkowski sum of two manifolds.
+ *
+ * @param other The other manifold to minkowski sum to this one.
+ * @param inset Whether it should add or subtract from the manifold.
+ */
+Manifold Manifold::Minkowski(const Manifold& other, bool inset) const {
+  std::vector<Manifold> composedHulls({*this});
+  bool aConvex = this->GetCsgLeafNode().GetImpl()->ReflexFaces().size() == 0;
+  bool bConvex = other.GetCsgLeafNode().GetImpl()->ReflexFaces().size() == 0;
+
+  // If the convex manifold was supplied first, swap them!
+  Manifold a = *this, b = other;
+  if (aConvex && !bConvex) {
+    a = other;
+    b = *this;
+    aConvex = !aConvex;
+    bConvex = !bConvex;
+  }
+
+  manifold::Mesh aMesh = a.GetMesh();
+
+  // Convex-Convex Minkowski: Very Fast
+  if (!inset && aConvex && bConvex) {
+    std::vector<Manifold> simpleHull;
+    for (glm::vec3 vertex : aMesh.vertPos) {
+      simpleHull.push_back(b.Translate(vertex));
+    }
+    composedHulls.push_back(Manifold::Hull(simpleHull));
+    // Convex - Non-Convex Minkowski: Slower
+  } else if ((inset || !aConvex) && bConvex) {
+    std::vector<std::vector<Manifold>> composedParts;
+    for (glm::ivec3 vertexIndices : aMesh.triVerts) {
+      composedParts.push_back({b.Translate(aMesh.vertPos[vertexIndices.x]),
+                               b.Translate(aMesh.vertPos[vertexIndices.y]),
+                               b.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 = b.GetMesh();
+    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);
+}
+
+/**
+ * Compute the minkowski sum of this manifold with another.
+ * This corresponds to the morphological dilation of the manifold.
+ *
+ * @param other The other manifold to minkowski sum to this one.
+ */
+Manifold Manifold::MinkowskiSum(const Manifold& other) const {
+  return this->Minkowski(other, false);
+}
+
+/**
+ * Subtract the sweep of the other manifold across this manifold's surface.
+ * This corresponds to the morphological erosion of the manifold.
+ *
+ * @param other The other manifold to minkowski subtract from this one.
+ */
+Manifold Manifold::MinkowskiDifference(const Manifold& other) const {
+  return this->Minkowski(other, true);
+}
+
 /**
  * Returns the cross section of this object parallel to the X-Y plane at the
  * specified Z height, defaulting to zero. Using a height equal to the bottom of

test/boolean_test.cpp

@@ -296,6 +296,67 @@ TEST(Boolean, SplitByPlane60) {
               splits.second.GetProperties().volume, 1e-5);
 }
 
+TEST(Manifold, ConvexConvexMinkowski) {
+  bool oldDeterministic = ManifoldParams().deterministic;
+  ManifoldParams().deterministic = true;
+  float offsetRadius = 0.1f;
+  float cubeWidth = 2.0f;
+  Manifold sphere = Manifold::Sphere(offsetRadius, 20);
+  Manifold cube = Manifold::Cube({cubeWidth, cubeWidth, cubeWidth});
+  Manifold sum = cube.MinkowskiSum(sphere);
+  EXPECT_NEAR(sum.GetProperties().volume,
+              powf(cubeWidth + (2 * offsetRadius), 3.f), 0.06f);
+  EXPECT_EQ(sum.Genus(), 0);
+  Manifold difference = Manifold::Cube({cubeWidth, cubeWidth, cubeWidth})
+                            .MinkowskiDifference(sphere);
+  EXPECT_NEAR(difference.GetProperties().volume,
+              powf(cubeWidth - (2 * offsetRadius), 3.f), 1e-5);
+  EXPECT_EQ(difference.Genus(), 0);
+  ManifoldParams().deterministic = oldDeterministic;
+}
+
+TEST(Manifold, NonConvexConvexMinkowski) {
+  bool oldDeterministic = ManifoldParams().deterministic;
+  ManifoldParams().deterministic = true;
+  Manifold sphere = Manifold::Sphere(1.2, 20);
+  Manifold cube = Manifold::Cube({2.0, 2.0, 2.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum = nonConvex.MinkowskiSum(Manifold::Sphere(0.1, 20));
+  EXPECT_NEAR(sum.GetProperties().volume, 4.8406339f, 1e-5);
+  EXPECT_EQ(sum.Genus(), 5);
+  Manifold difference =
+      nonConvex.MinkowskiDifference(Manifold::Sphere(0.05, 20));
+  EXPECT_NEAR(difference.GetProperties().volume, 0.77841246128082275f, 1e-5);
+  EXPECT_EQ(difference.Genus(), 5);
+  ManifoldParams().deterministic = oldDeterministic;
+}
+
+TEST(Manifold, DISABLED_NonConvexNonConvexMinkowski) {
+  bool oldDeterministic = ManifoldParams().deterministic;
+  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(1.2, 20);
+  Manifold cube = Manifold::Cube({2.0, 2.0, 2.0}, true);
+  Manifold nonConvex = cube - sphere;
+  Manifold sum = nonConvex.MinkowskiSum(star);
+  EXPECT_NEAR(sum.GetProperties().volume, 7.0875549f, 1e-5);
+  EXPECT_EQ(sum.Genus(), 5);
+  Manifold difference = nonConvex.MinkowskiDifference(star);
+  EXPECT_NEAR(difference.GetProperties().volume, 0.0093147634f, 1e-5);
+  EXPECT_EQ(difference.Genus(), -7);
+  ManifoldParams().deterministic = oldDeterministic;
+}
+
 /**
  * This tests that non-intersecting geometry is properly retained.
  */