diff --git a/.vscode/settings.json b/.vscode/settings.json index 8fc9dbdf2..3b43ef926 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -146,6 +146,7 @@ "Gyroid", "halfedge", "halfedges", + "Minkowski", "Tris", "Verts", "Voxel", diff --git a/bindings/c/include/manifold/manifoldc.h b/bindings/c/include/manifold/manifoldc.h index fdb330d00..9bfa856b8 100644 --- a/bindings/c/include/manifold/manifoldc.h +++ b/bindings/c/include/manifold/manifoldc.h @@ -115,6 +115,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 diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp index 78795eb1a..d862f3d23 100644 --- a/bindings/c/manifoldc.cpp +++ b/bindings/c/manifoldc.cpp @@ -183,6 +183,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); diff --git a/bindings/python/examples/minkowski.py b/bindings/python/examples/minkowski.py new file mode 100644 index 000000000..2cbbda471 --- /dev/null +++ b/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) diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index cb5e6f8b8..8f2e1edb2 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -366,6 +366,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) { diff --git a/bindings/wasm/bindings.cpp b/bindings/wasm/bindings.cpp index c396b4e42..77cb9d2fe 100644 --- a/bindings/wasm/bindings.cpp +++ b/bindings/wasm/bindings.cpp @@ -155,6 +155,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::Hull)) diff --git a/bindings/wasm/examples/worker.ts b/bindings/wasm/examples/worker.ts index ace0e8524..ffee1e45f 100644 --- a/bindings/wasm/examples/worker.ts +++ b/bindings/wasm/examples/worker.ts @@ -83,7 +83,9 @@ const manifoldMemberFunctions = [ 'splitByPlane', 'slice', 'project', - 'hull' + 'hull', + 'minkowskiSum', + 'minkowskiDifference' ]; // CrossSection static methods (that return a new cross-section) const crossSectionStaticFunctions = [ diff --git a/bindings/wasm/manifold-encapsulated-types.d.ts b/bindings/wasm/manifold-encapsulated-types.d.ts index f9154d2af..f2d3b8775 100644 --- a/bindings/wasm/manifold-encapsulated-types.d.ts +++ b/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 diff --git a/include/manifold/manifold.h b/include/manifold/manifold.h index 9c59d7367..590584733 100644 --- a/include/manifold/manifold.h +++ b/include/manifold/manifold.h @@ -331,6 +331,8 @@ class Manifold { std::pair SplitByPlane(vec3 normal, double originOffset) const; Manifold TrimByPlane(vec3 normal, double originOffset) const; + Manifold MinkowskiSum(const Manifold&) const; + Manifold MinkowskiDifference(const Manifold&) const; ///@} /** @name Properties @@ -391,6 +393,8 @@ class Manifold { mutable std::shared_ptr pNode_; CsgLeafNode& GetCsgLeafNode() const; + + Manifold Minkowski(const Manifold&, bool inset = false) const; }; /** @} */ diff --git a/src/edge_op.cpp b/src/edge_op.cpp index d22577a2a..44c8dcbbb 100644 --- a/src/edge_op.cpp +++ b/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_[idx / 3]; + const vec3 normal1 = pImpl->faceNormal_[edge.pairedHalfedge / 3]; + + if (linalg::all(linalg::equal(normal0, normal1))) return; + + const vec3 edgeVec = + pImpl->vertPos_[edge.endVert] - pImpl->vertPos_[edge.startVert]; + const bool convex = + linalg::dot(edgeVec, linalg::cross(normal0, normal1)) > tolerance; + if (!convex) anyConcave = true; + }); + + return !anyConcave; +} } // namespace manifold diff --git a/src/impl.h b/src/impl.h index d37a9a7d2..d464aa07a 100644 --- a/src/impl.h +++ b/src/impl.h @@ -320,6 +320,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; diff --git a/src/manifold.cpp b/src/manifold.cpp index a36dd460b..a4e707743 100644 --- a/src/manifold.cpp +++ b/src/manifold.cpp @@ -968,6 +968,107 @@ 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 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 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 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 = + linalg::all(linalg::equal(aImpl->faceNormal_[aFace], + bImpl->faceNormal_[bFace])) || + linalg::all(linalg::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 diff --git a/test/boolean_test.cpp b/test/boolean_test.cpp index 40bb00dd3..a931984a6 100644 --- a/test/boolean_test.cpp +++ b/test/boolean_test.cpp @@ -344,6 +344,75 @@ TEST(Boolean, SplitByPlane60) { EXPECT_NEAR(splits.first.Volume(), splits.second.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.Volume(), 10.589364051818848f, 1e-5); + EXPECT_EQ(sum.Genus(), 0); + Manifold difference = Manifold::Cube({cubeWidth, cubeWidth, cubeWidth}) + .MinkowskiDifference(sphere); + EXPECT_NEAR(difference.Volume(), 5.8319993019104004f, 1e-5); + EXPECT_NEAR(difference.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) { + 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.Volume(), 4.8406339f, 1e-5); + EXPECT_NEAR(sum.SurfaceArea(), 34.063014984130859f, 1e-5); + EXPECT_EQ(sum.Genus(), 5); + Manifold difference = + nonConvex.MinkowskiDifference(Manifold::Sphere(0.05, 20)); + EXPECT_NEAR(difference.Volume(), 0.77841246128082275f, 1e-5); + EXPECT_NEAR(difference.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().processOverlaps = false; +} + +TEST(Boolean, NonConvexNonConvexMinkowski) { + 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.Volume(), 8.65625f, 1e-5); + EXPECT_NEAR(sum.SurfaceArea(), 31.176914f, 1e-5); + EXPECT_EQ(sum.Genus(), 0); + + Manifold difference = + nonConvex.MinkowskiDifference(nonConvex.Scale(vec3(0.1))); + EXPECT_NEAR(difference.Volume(), 0.81554f, 1e-5); + EXPECT_NEAR(difference.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().processOverlaps = false; +} + /** * This tests that non-intersecting geometry is properly retained. */