diff --git a/.vscode/settings.json b/.vscode/settings.json index c229af2f8..8eb2bc012 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -147,6 +147,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 79597157c..fcc3faa00 100644 --- a/bindings/c/include/manifold/manifoldc.h +++ b/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 diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp index 1569bc061..ab649f41a 100644 --- a/bindings/c/manifoldc.cpp +++ b/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); 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 949b66ede..d13bea3e0 100644 --- a/bindings/python/manifold3d.cpp +++ b/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) { diff --git a/bindings/wasm/bindings.cpp b/bindings/wasm/bindings.cpp index 881ace443..25cfb28f6 100644 --- a/bindings/wasm/bindings.cpp +++ b/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::Hull)) diff --git a/bindings/wasm/examples/worker.ts b/bindings/wasm/examples/worker.ts index d8151d9b1..8443d1d45 100644 --- a/bindings/wasm/examples/worker.ts +++ b/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 = [ diff --git a/bindings/wasm/manifold-encapsulated-types.d.ts b/bindings/wasm/manifold-encapsulated-types.d.ts index 24e9ce4fb..0e2ca4521 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 f3007703c..64d0a5c31 100644 --- a/include/manifold/manifold.h +++ b/include/manifold/manifold.h @@ -290,6 +290,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 2D from 3D @@ -325,6 +327,8 @@ class Manifold { mutable std::shared_ptr pNode_; CsgLeafNode& GetCsgLeafNode() const; + + Manifold Minkowski(const Manifold&, bool inset = false) const; }; /** @} */ } // namespace manifold diff --git a/src/edge_op.cpp b/src/edge_op.cpp index 759a81031..2b912e533 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_[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 diff --git a/src/impl.h b/src/impl.h index 45e3c95f4..d47083f39 100644 --- a/src/impl.h +++ b/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; diff --git a/src/manifold.cpp b/src/manifold.cpp index f29387950..a8c0baba9 100644 --- a/src/manifold.cpp +++ b/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 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 = 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 diff --git a/test/boolean_test.cpp b/test/boolean_test.cpp index 043768c40..b660778cf 100644 --- a/test/boolean_test.cpp +++ b/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. */