Add MinkowskiSum and MinkowskiDifference

.vscode/settings.json

@@ -147,6 +147,7 @@
     "Gyroid",
     "halfedge",
     "halfedges",
+    "Minkowski",
     "Tris",
     "Verts",
     "Voxel",

bindings/c/include/manifold/manifoldc.h

@@ -104,6 +104,10 @@ ManifoldManifoldPair manifold_split_by_plane(void *mem_first, void *mem_second,
 ManifoldManifold *manifold_trim_by_plane(void *mem, ManifoldManifold *m,
                                          double normal_x, double normal_y,
                                          double normal_z, double 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

@@ -182,6 +182,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));
+}
+
 ManifoldPolygons *manifold_slice(void *mem, ManifoldManifold *m,
                                  double 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, 15)
+    cube = Manifold.cube([1.0, 1.0, 1.0], True)
+    sphereless_cube = cube - sphere
+    return sphereless_cube.minkowski_sum(star)

bindings/python/manifold3d.cpp

@@ -363,6 +363,10 @@ NB_MODULE(manifold3d, m) {
       .def("trim_by_plane", &Manifold::TrimByPlane, nb::arg("normal"),
            nb::arg("origin_offset"),
            manifold__trim_by_plane__normal__origin_offset)
+      .def("minkowski_sum", &Manifold::MinkowskiSum, nb::arg("other"),
+           manifold__minkowski_sum__other)
+      .def("minkowski_difference", &Manifold::MinkowskiDifference,
+           nb::arg("other"), manifold__minkowski_difference__other)
       .def(
           "slice",
           [](const Manifold &self, double height) {

bindings/wasm/bindings.cpp

@@ -135,6 +135,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

@@ -82,7 +82,9 @@ const manifoldMemberFunctions = [
   'splitByPlane',
   'slice',
   'project',
-  'hull'
+  'hull',
+  'minkowskiSum',
+  'minkowskiDifference'
 ];
 // CrossSection static methods (that return a new cross-section)
 const crossSectionStaticFunctions = [

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

@@ -841,6 +841,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

include/manifold/manifold.h

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

src/edge_op.cpp

@@ -693,4 +693,36 @@ void Manifold::Impl::SplitPinchedVerts() {
     });
   }
 }
+
+// Return true if Manifold is Genus 0 and contains no concave edges
+bool Manifold::Impl::IsConvex(float tolerance) const {
+  // Convex Shape Must have Genus of 0
+  int chi = NumVert() - NumEdge() + NumTri();
+  int genus = 1 - chi / 2;
+  if (genus != 0) return false;
+
+  // Iterate across all edges; return false if any edges are concave
+  const Impl* pImpl = this;
+  const size_t nbEdges = halfedge_.size();
+  auto policy = autoPolicy(nbEdges, 1e5);
+  bool anyConcave = false;
+  for_each_n(
+      policy, countAt(0), nbEdges, [&anyConcave, &pImpl, &tolerance](int idx) {
+        Halfedge edge = pImpl->halfedge_[idx];
+        if (!edge.IsForward()) return;
+
+        const vec3 normal0 = pImpl->faceNormal_[edge.face];
+        const vec3 normal1 =
+            pImpl->faceNormal_[pImpl->halfedge_[edge.pairedHalfedge].face];
+        if (glm::all(glm::equal(normal0, normal1))) return;
+
+        const vec3 edgeVec =
+            pImpl->vertPos_[edge.endVert] - pImpl->vertPos_[edge.startVert];
+        const bool convex =
+            glm::dot(edgeVec, glm::cross(normal0, normal1)) > tolerance;
+        if (!convex) anyConcave = true;
+      });
+
+  return !anyConcave;
+}
 }  // namespace manifold

src/impl.h

@@ -310,6 +310,7 @@ struct Manifold::Impl {
   void FormLoop(int current, int end);
   void CollapseTri(const ivec3& triEdge);
   void SplitPinchedVerts();
+  bool IsConvex(float tolerance = 1e-8f) const;
 
   // subdivision.cpp
   int GetNeighbor(int tri) const;

src/manifold.cpp

@@ -945,6 +945,106 @@ Manifold Manifold::TrimByPlane(vec3 normal, double originOffset) const {
   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});
+  auto aImpl = this->GetCsgLeafNode().GetImpl();
+  auto bImpl = other.GetCsgLeafNode().GetImpl();
+
+  bool aConvex = aImpl->IsConvex();
+  bool bConvex = bImpl->IsConvex();
+
+  // If the convex manifold was supplied first, swap them!
+  Manifold a = *this, b = other;
+  if (aConvex && !bConvex) {
+    a = other;
+    b = *this;
+    aImpl = other.GetCsgLeafNode().GetImpl();
+    bImpl = this->GetCsgLeafNode().GetImpl();
+    aConvex = !aConvex;
+    bConvex = !bConvex;
+  }
+
+  // Early-exit if either input is empty
+  if (b.IsEmpty()) return a;
+  if (a.IsEmpty()) return b;
+
+  // Convex-Convex Minkowski: Very Fast
+  if (!inset && aConvex && bConvex) {
+    std::vector<Manifold> simpleHull;
+    for (const vec3& vertex : aImpl->vertPos_) {
+      simpleHull.push_back(b.Translate(vertex));
+    }
+    composedHulls.push_back(Manifold::Hull(simpleHull));
+    // Convex - Non-Convex Minkowski: Slower
+  } else if ((inset || !aConvex) && bConvex) {
+    const size_t numTri = aImpl->NumTri();
+    std::vector<Manifold> newHulls(numTri);
+    auto policy = autoPolicy(numTri, 100);
+    for_each_n(
+        policy, countAt(0), numTri, [&newHulls, &b, &aImpl](const int face) {
+          newHulls[face] = Manifold::Hull(
+              {b.Translate(
+                   aImpl->vertPos_[aImpl->halfedge_[(face * 3) + 0].startVert]),
+               b.Translate(
+                   aImpl->vertPos_[aImpl->halfedge_[(face * 3) + 1].startVert]),
+               b.Translate(aImpl->vertPos_[aImpl->halfedge_[(face * 3) + 2]
+                                               .startVert])});
+        });
+    composedHulls.insert(composedHulls.end(), newHulls.begin(), newHulls.end());
+    // Non-Convex - Non-Convex Minkowski: Very Slow
+  } else if (!aConvex && !bConvex) {
+    for (size_t aFace = 0; aFace < aImpl->NumTri(); aFace++) {
+      for (size_t bFace = 0; bFace < bImpl->NumTri(); bFace++) {
+        const bool coplanar = glm::all(glm::equal(aImpl->faceNormal_[aFace],
+                                                  bImpl->faceNormal_[bFace])) ||
+                              glm::all(glm::equal(aImpl->faceNormal_[aFace],
+                                                  -bImpl->faceNormal_[bFace]));
+        if (coplanar) continue;  // Skip Coplanar Triangles
+
+        vec3 a1 = aImpl->vertPos_[aImpl->halfedge_[(aFace * 3) + 0].startVert];
+        vec3 a2 = aImpl->vertPos_[aImpl->halfedge_[(aFace * 3) + 1].startVert];
+        vec3 a3 = aImpl->vertPos_[aImpl->halfedge_[(aFace * 3) + 2].startVert];
+        vec3 b1 = bImpl->vertPos_[bImpl->halfedge_[(bFace * 3) + 0].startVert];
+        vec3 b2 = bImpl->vertPos_[bImpl->halfedge_[(bFace * 3) + 1].startVert];
+        vec3 b3 = bImpl->vertPos_[bImpl->halfedge_[(bFace * 3) + 2].startVert];
+        composedHulls.push_back(
+            Manifold::Hull({a1 + b1, a1 + b2, a1 + b3, a2 + b1, a2 + b2,
+                            a2 + b3, a3 + b1, a3 + b2, a3 + b3}));
+      }
+    }
+  }
+  return Manifold::BatchBoolean(composedHulls, inset
+                                                   ? manifold::OpType::Subtract
+                                                   : manifold::OpType::Add)
+      .AsOriginal();
+}
+
+/**
+ * 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

@@ -379,6 +379,81 @@ TEST(Boolean, SplitByPlane60) {
               splits.second.GetProperties().volume, 1e-5);
 }
 
+TEST(Boolean, ConvexConvexMinkowski) {
+  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, 10.589364051818848f, 1e-5);
+  EXPECT_EQ(sum.Genus(), 0);
+  Manifold difference = Manifold::Cube({cubeWidth, cubeWidth, cubeWidth})
+                            .MinkowskiDifference(sphere);
+  EXPECT_NEAR(difference.GetProperties().volume, 5.8319993019104004f, 1e-5);
+  EXPECT_NEAR(difference.GetProperties().surfaceArea, 19.439998626708984, 1e-5);
+  EXPECT_EQ(difference.Genus(), 0);
+
+#ifdef MANIFOLD_EXPORT
+  if (options.exportModels)
+    ExportMesh("minkowski-convex-convex.glb", sum.GetMeshGL(), {});
+#endif
+}
+
+TEST(Boolean, NonConvexConvexMinkowski) {
+  bool oldDeterministic = ManifoldParams().deterministic;
+  ManifoldParams().deterministic = true;
+  ManifoldParams().processOverlaps = 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_NEAR(sum.GetProperties().surfaceArea, 34.063014984130859f, 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_NEAR(difference.GetProperties().surfaceArea, 16.703740785913258, 1e-5);
+  EXPECT_EQ(difference.Genus(), 5);
+
+#ifdef MANIFOLD_EXPORT
+  if (options.exportModels)
+    ExportMesh("minkowski-nonconvex-convex.glb", sum.GetMeshGL(), {});
+#endif
+
+  ManifoldParams().deterministic = oldDeterministic;
+  ManifoldParams().processOverlaps = false;
+}
+
+TEST(Boolean, NonConvexNonConvexMinkowski) {
+  bool oldDeterministic = ManifoldParams().deterministic;
+  ManifoldParams().deterministic = true;
+  ManifoldParams().processOverlaps = true;
+
+  Manifold tet = Manifold::Tetrahedron();
+  Manifold nonConvex = tet - tet.Rotate(0, 0, 90).Translate(vec3(1));
+
+  Manifold sum = nonConvex.MinkowskiSum(nonConvex.Scale(vec3(0.5)));
+  EXPECT_NEAR(sum.GetProperties().volume, 8.65625f, 1e-5);
+  EXPECT_NEAR(sum.GetProperties().surfaceArea, 31.176914f, 1e-5);
+  EXPECT_EQ(sum.Genus(), 0);
+
+  Manifold difference =
+      nonConvex.MinkowskiDifference(nonConvex.Scale(vec3(0.1)));
+  EXPECT_NEAR(difference.GetProperties().volume, 0.81554f, 1e-5);
+  EXPECT_NEAR(difference.GetProperties().surfaceArea, 6.95045f, 1e-5);
+  EXPECT_EQ(difference.Genus(), 0);
+
+#ifdef MANIFOLD_EXPORT
+  if (options.exportModels)
+    ExportMesh("minkowski-nonconvex-nonconvex.glb", sum.GetMeshGL(), {});
+#endif
+
+  ManifoldParams().deterministic = oldDeterministic;
+  ManifoldParams().processOverlaps = false;
+}
+
 /**
  * This tests that non-intersecting geometry is properly retained.
  */