diff --git a/.vscode/settings.json b/.vscode/settings.json index 378812de1..b3ae887d5 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -150,4 +150,4 @@ }, "python.formatting.provider": "none", "clang-format.executable": "clang-format", -} \ No newline at end of file +} diff --git a/bindings/c/include/manifoldc.h b/bindings/c/include/manifoldc.h index 18cc408b4..88a8161f8 100644 --- a/bindings/c/include/manifoldc.h +++ b/bindings/c/include/manifoldc.h @@ -190,6 +190,8 @@ ManifoldManifold *manifold_set_properties( void (*fun)(float *new_prop, ManifoldVec3 position, const float *old_prop)); ManifoldManifold *manifold_calculate_curvature(void *mem, ManifoldManifold *m, int gaussian_idx, int mean_idx); +float manifold_min_gap(ManifoldManifold *m, ManifoldManifold *other, + float searchLength); ManifoldManifold *manifold_calculate_normals(void *mem, ManifoldManifold *m, int normal_idx, int min_sharp_angle); diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp index 895b3f127..521557069 100644 --- a/bindings/c/manifoldc.cpp +++ b/bindings/c/manifoldc.cpp @@ -550,6 +550,11 @@ ManifoldManifold *manifold_calculate_curvature(void *mem, ManifoldManifold *m, return to_c(new (mem) Manifold(man)); } +float manifold_min_gap(ManifoldManifold *m, ManifoldManifold *other, + float searchLength) { + return from_c(m)->MinGap(*from_c(other), searchLength); +} + ManifoldManifold *manifold_calculate_normals(void *mem, ManifoldManifold *m, int normal_idx, int min_sharp_angle) { diff --git a/bindings/python/examples/all_apis.py b/bindings/python/examples/all_apis.py index 53e52fab0..c2e231fb5 100644 --- a/bindings/python/examples/all_apis.py +++ b/bindings/python/examples/all_apis.py @@ -111,6 +111,9 @@ def all_manifold(): m = m.trim_by_plane((0, 0, 1), 0) m = m.warp(lambda p: (p[0] + 1, p[1] / 2, p[2] * 2)) m = m.warp_batch(lambda ps: ps * [1, 0.5, 2] + [1, 0, 0]) + m = Manifold.cube() + m2 = Manifold.cube().translate([2, 0, 0]) + d = m.min_gap(m2, 2) def run(): diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index 768a4cd9f..c9f3a5efd 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -315,6 +315,10 @@ NB_MODULE(manifold3d, m) { .def("calculate_curvature", &Manifold::CalculateCurvature, nb::arg("gaussian_idx"), nb::arg("mean_idx"), manifold__calculate_curvature__gaussian_idx__mean_idx) + .def("min_gap", &Manifold::MinGap, nb::arg("other"), + nb::arg("search_length"), + "Returns the minimum gap between two manifolds." + "Returns a float between 0 and searchLength.") .def("calculate_normals", &Manifold::CalculateNormals, nb::arg("normal_idx"), nb::arg("min_sharp_angle") = 60, manifold__calculate_normals__normal_idx__min_sharp_angle) diff --git a/bindings/wasm/bindings.cpp b/bindings/wasm/bindings.cpp index 5c703e531..de1e7b22a 100644 --- a/bindings/wasm/bindings.cpp +++ b/bindings/wasm/bindings.cpp @@ -165,6 +165,7 @@ EMSCRIPTEN_BINDINGS(whatever) { .function("precision", &Manifold::Precision) .function("genus", &Manifold::Genus) .function("getProperties", &Manifold::GetProperties) + .function("minGap", &Manifold::MinGap) .function("calculateCurvature", &Manifold::CalculateCurvature) .function("calculateNormals", &Manifold::CalculateNormals) .function("originalID", &Manifold::OriginalID) diff --git a/bindings/wasm/manifold-encapsulated-types.d.ts b/bindings/wasm/manifold-encapsulated-types.d.ts index 8482d3fdc..d51afcae9 100644 --- a/bindings/wasm/manifold-encapsulated-types.d.ts +++ b/bindings/wasm/manifold-encapsulated-types.d.ts @@ -845,6 +845,13 @@ export class Manifold { */ getProperties(): Properties; + + /* + * Returns the minimum gap between two manifolds. Returns a float between + * 0 and searchLength. + */ + minGap(other: Manifold, searchLength: number): number; + // Export /** diff --git a/src/manifold/include/manifold.h b/src/manifold/include/manifold.h index a99006c4b..a79cbbc7b 100644 --- a/src/manifold/include/manifold.h +++ b/src/manifold/include/manifold.h @@ -218,6 +218,7 @@ class Manifold { float Precision() const; int Genus() const; Properties GetProperties() const; + float MinGap(const Manifold& other, float searchLength) const; ///@} /** @name Mesh ID diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index 8de196a08..014d2945b 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -23,6 +23,7 @@ #include "mesh_fixes.h" #include "par.h" #include "svd.h" +#include "tri_dist.h" namespace { using namespace manifold; @@ -902,4 +903,48 @@ SparseIndices Manifold::Impl::VertexCollisionsZ( else return collider_.Collisions(vertsIn); } + +/* + * Returns the minimum gap between two manifolds. Returns a float between + * 0 and searchLength. + */ +float Manifold::Impl::MinGap(const Manifold::Impl& other, + float searchLength) const { + ZoneScoped; + Vec faceBoxOther; + Vec faceMortonOther; + + other.GetFaceBoxMorton(faceBoxOther, faceMortonOther); + + transform(autoPolicy(faceBoxOther.size()), faceBoxOther.begin(), + faceBoxOther.end(), faceBoxOther.begin(), + [searchLength](const Box& box) { + return Box(box.min - glm::vec3(searchLength), + box.max + glm::vec3(searchLength)); + }); + + SparseIndices collisions = collider_.Collisions(faceBoxOther.cview()); + + float minDistanceSquared = transform_reduce( + autoPolicy(collisions.size()), thrust::counting_iterator(0), + thrust::counting_iterator(collisions.size()), + [&collisions, this, &other](int i) { + const int tri = collisions.Get(i, 1); + const int triOther = collisions.Get(i, 0); + + std::array p; + std::array q; + + for (const int j : {0, 1, 2}) { + p[j] = vertPos_[halfedge_[3 * tri + j].startVert]; + q[j] = other.vertPos_[other.halfedge_[3 * triOther + j].startVert]; + } + + return DistanceTriangleTriangleSquared(p, q); + }, + searchLength * searchLength, thrust::minimum()); + + return sqrt(minDistanceSquared); +}; + } // namespace manifold diff --git a/src/manifold/src/impl.h b/src/manifold/src/impl.h index d0f94a8a2..6f026b730 100644 --- a/src/manifold/src/impl.h +++ b/src/manifold/src/impl.h @@ -81,6 +81,7 @@ struct Manifold::Impl { SparseIndices EdgeCollisions(const Impl& B, bool inverted = false) const; SparseIndices VertexCollisionsZ(VecView vertsIn, bool inverted = false) const; + float MinGap(const Impl& other, float searchLength) const; bool IsEmpty() const { return NumVert() == 0; } int NumVert() const { return vertPos_.size(); } diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 797a42b08..810405d46 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -21,6 +21,7 @@ #include "csg_tree.h" #include "impl.h" #include "par.h" +#include "tri_dist.h" namespace { using namespace manifold; @@ -926,4 +927,21 @@ Manifold Manifold::Hull() const { return Hull(GetMesh().vertPos); } Manifold Manifold::Hull(const std::vector& manifolds) { return Compose(manifolds).Hull(); } + +/** + * Returns the minimum gap between two manifolds. Returns a float between + * 0 and searchLength. + * + * @param other The other manifold to compute the minimum gap to. + * @param searchLength The maximum distance to search for a minimum gap. + */ +float Manifold::MinGap(const Manifold& other, float searchLength) const { + auto intersect = *this ^ other; + auto prop = intersect.GetProperties(); + + if (prop.volume != 0) return 0.0f; + + return GetCsgLeafNode().GetImpl()->MinGap(*other.GetCsgLeafNode().GetImpl(), + searchLength); +} } // namespace manifold diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index bdd2253a3..83e7b6ca2 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -83,4 +83,4 @@ endif() target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_17) install(TARGETS ${PROJECT_NAME} EXPORT manifoldTargets) -install(FILES include/public.h include/vec_view.h DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}) +install(FILES include/public.h include/vec_view.h include/tri_dist.h DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}) diff --git a/src/utilities/include/tri_dist.h b/src/utilities/include/tri_dist.h new file mode 100644 index 000000000..f768f7f90 --- /dev/null +++ b/src/utilities/include/tri_dist.h @@ -0,0 +1,225 @@ +// Copyright 2024 The Manifold Authors. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include + +namespace manifold { + +// From NVIDIA-Omniverse PhysX - BSD 3-Clause "New" or "Revised" License +// https://github.com/NVIDIA-Omniverse/PhysX/blob/main/LICENSE.md +// https://github.com/NVIDIA-Omniverse/PhysX/blob/main/physx/source/geomutils/src/sweep/GuSweepCapsuleCapsule.cpp +// With minor modifications + +/** + * Returns the distance between two line segments. + * + * @param[out] x Closest point on line segment pa. + * @param[out] y Closest point on line segment qb. + * @param[in] p One endpoint of the first line segment. + * @param[in] a Other endpoint of the first line segment. + * @param[in] p One endpoint of the second line segment. + * @param[in] b Other endpoint of the second line segment. + */ +inline void EdgeEdgeDist(glm::vec3& x, glm::vec3& y, // closest points + const glm::vec3& p, + const glm::vec3& a, // seg 1 origin, vector + const glm::vec3& q, + const glm::vec3& b) // seg 2 origin, vector +{ + const glm::vec3 T = q - p; + const float ADotA = glm::dot(a, a); + const float BDotB = glm::dot(b, b); + const float ADotB = glm::dot(a, b); + const float ADotT = glm::dot(a, T); + const float BDotT = glm::dot(b, T); + + // t parameterizes ray (p, a) + // u parameterizes ray (q, b) + + // Compute t for the closest point on ray (p, a) to ray (q, b) + const float Denom = ADotA * BDotB - ADotB * ADotB; + + float t; // We will clamp result so t is on the segment (p, a) + t = Denom != 0.0f + ? glm::clamp((ADotT * BDotB - BDotT * ADotB) / Denom, 0.0f, 1.0f) + : 0.0f; + + // find u for point on ray (q, b) closest to point at t + float u; + if (BDotB != 0.0f) { + u = (t * ADotB - BDotT) / BDotB; + + // if u is on segment (q, b), t and u correspond to closest points, + // otherwise, clamp u, recompute and clamp t + if (u < 0.0f) { + u = 0.0f; + t = ADotA != 0.0f ? glm::clamp(ADotT / ADotA, 0.0f, 1.0f) : 0.0f; + } else if (u > 1.0f) { + u = 1.0f; + t = ADotA != 0.0f ? glm::clamp((ADotB + ADotT) / ADotA, 0.0f, 1.0f) + : 0.0f; + } + } else { + u = 0.0f; + t = ADotA != 0.0f ? glm::clamp(ADotT / ADotA, 0.0f, 1.0f) : 0.0f; + } + x = p + a * t; + y = q + b * u; +} + +// From NVIDIA-Omniverse PhysX - BSD 3-Clause "New" or "Revised" License +// https://github.com/NVIDIA-Omniverse/PhysX/blob/main/LICENSE.md +// https://github.com/NVIDIA-Omniverse/PhysX/blob/main/physx/source/geomutils/src/distance/GuDistanceTriangleTriangle.cpp +// With minor modifications + +/** + * Returns the minimum squared distance between two triangles. + * + * @param p First triangle. + * @param q Second triangle. + */ +inline float DistanceTriangleTriangleSquared( + const std::array& p, const std::array& q) { + std::array Sv; + Sv[0] = p[1] - p[0]; + Sv[1] = p[2] - p[1]; + Sv[2] = p[0] - p[2]; + + std::array Tv; + Tv[0] = q[1] - q[0]; + Tv[1] = q[2] - q[1]; + Tv[2] = q[0] - q[2]; + + bool shown_disjoint = false; + + float mindd = std::numeric_limits::max(); + + for (uint32_t i = 0; i < 3; i++) { + for (uint32_t j = 0; j < 3; j++) { + glm::vec3 cp; + glm::vec3 cq; + EdgeEdgeDist(cp, cq, p[i], Sv[i], q[j], Tv[j]); + const glm::vec3 V = cq - cp; + const float dd = glm::dot(V, V); + + if (dd <= mindd) { + mindd = dd; + + uint32_t id = i + 2; + if (id >= 3) id -= 3; + glm::vec3 Z = p[id] - cp; + float a = glm::dot(Z, V); + id = j + 2; + if (id >= 3) id -= 3; + Z = q[id] - cq; + float b = glm::dot(Z, V); + + if ((a <= 0.0f) && (b >= 0.0f)) { + return glm::dot(V, V); + }; + + if (a <= 0.0f) + a = 0.0f; + else if (b > 0.0f) + b = 0.0f; + + if ((mindd - a + b) > 0.0f) shown_disjoint = true; + } + } + } + + glm::vec3 Sn = glm::cross(Sv[0], Sv[1]); + float Snl = glm::dot(Sn, Sn); + + if (Snl > 1e-15f) { + const glm::vec3 Tp(glm::dot(p[0] - q[0], Sn), glm::dot(p[0] - q[1], Sn), + glm::dot(p[0] - q[2], Sn)); + + int index = -1; + if ((Tp[0] > 0.0f) && (Tp[1] > 0.0f) && (Tp[2] > 0.0f)) { + index = Tp[0] < Tp[1] ? 0 : 1; + if (Tp[2] < Tp[index]) index = 2; + } else if ((Tp[0] < 0.0f) && (Tp[1] < 0.0f) && (Tp[2] < 0.0f)) { + index = Tp[0] > Tp[1] ? 0 : 1; + if (Tp[2] > Tp[index]) index = 2; + } + + if (index >= 0) { + shown_disjoint = true; + + const glm::vec3& qIndex = q[index]; + + glm::vec3 V = qIndex - p[0]; + glm::vec3 Z = glm::cross(Sn, Sv[0]); + if (glm::dot(V, Z) > 0.0f) { + V = qIndex - p[1]; + Z = glm::cross(Sn, Sv[1]); + if (glm::dot(V, Z) > 0.0f) { + V = qIndex - p[2]; + Z = glm::cross(Sn, Sv[2]); + if (glm::dot(V, Z) > 0.0f) { + glm::vec3 cp = qIndex + Sn * Tp[index] / Snl; + glm::vec3 cq = qIndex; + return glm::dot(cp - cq, cp - cq); + } + } + } + } + } + + glm::vec3 Tn = glm::cross(Tv[0], Tv[1]); + float Tnl = glm::dot(Tn, Tn); + + if (Tnl > 1e-15f) { + const glm::vec3 Sp(glm::dot(q[0] - p[0], Tn), glm::dot(q[0] - p[1], Tn), + glm::dot(q[0] - p[2], Tn)); + + int index = -1; + if ((Sp[0] > 0.0f) && (Sp[1] > 0.0f) && (Sp[2] > 0.0f)) { + index = Sp[0] < Sp[1] ? 0 : 1; + if (Sp[2] < Sp[index]) index = 2; + } else if ((Sp[0] < 0.0f) && (Sp[1] < 0.0f) && (Sp[2] < 0.0f)) { + index = Sp[0] > Sp[1] ? 0 : 1; + if (Sp[2] > Sp[index]) index = 2; + } + + if (index >= 0) { + shown_disjoint = true; + + const glm::vec3& pIndex = p[index]; + + glm::vec3 V = pIndex - q[0]; + glm::vec3 Z = glm::cross(Tn, Tv[0]); + if (glm::dot(V, Z) > 0.0f) { + V = pIndex - q[1]; + Z = glm::cross(Tn, Tv[1]); + if (glm::dot(V, Z) > 0.0f) { + V = pIndex - q[2]; + Z = glm::cross(Tn, Tv[2]); + if (glm::dot(V, Z) > 0.0f) { + glm::vec3 cp = pIndex; + glm::vec3 cq = pIndex + Tn * Sp[index] / Tnl; + return glm::dot(cp - cq, cp - cq); + } + } + } + } + } + + return shown_disjoint ? mindd : 0.0f; +}; +} // namespace manifold diff --git a/test/manifold_test.cpp b/test/manifold_test.cpp index f67bee9f3..5160c6283 100644 --- a/test/manifold_test.cpp +++ b/test/manifold_test.cpp @@ -17,7 +17,9 @@ #include #include "cross_section.h" +#include "samples.h" #include "test.h" +#include "tri_dist.h" #ifdef MANIFOLD_EXPORT #include "meshIO.h" @@ -833,3 +835,157 @@ TEST(Manifold, EmptyHull) { {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}}; EXPECT_TRUE(Manifold::Hull(coplanar).IsEmpty()); } + +TEST(Manifold, MinGapCubeCube) { + auto a = Manifold::Cube(); + auto b = Manifold::Cube().Translate({2, 2, 0}); + + float distance = a.MinGap(b, 1.5f); + + EXPECT_FLOAT_EQ(distance, sqrt(2)); +} + +TEST(Manifold, MinGapCubeCube2) { + auto a = Manifold::Cube(); + auto b = Manifold::Cube().Translate({3, 3, 0}); + + float distance = a.MinGap(b, 3); + + EXPECT_FLOAT_EQ(distance, sqrt(2) * 2); +} + +TEST(Manifold, MinGapCubeSphereOverlapping) { + auto a = Manifold::Cube(); + auto b = Manifold::Sphere(1); + + float distance = a.MinGap(b, 0.1f); + + EXPECT_FLOAT_EQ(distance, 0); +} + +TEST(Manifold, MinGapSphereSphere) { + auto a = Manifold::Sphere(1); + auto b = Manifold::Sphere(1).Translate({2, 2, 0}); + + float distance = a.MinGap(b, 0.85f); + + EXPECT_FLOAT_EQ(distance, 2 * sqrt(2) - 2); +} + +TEST(Manifold, MinGapSphereSphereOutOfBounds) { + auto a = Manifold::Sphere(1); + auto b = Manifold::Sphere(1).Translate({2, 2, 0}); + + float distance = a.MinGap(b, 0.8f); + + EXPECT_FLOAT_EQ(distance, 0.8f); +} + +TEST(Manifold, MinGapClosestPointOnEdge) { + auto a = Manifold::Cube({1, 1, 1}, true).Rotate(0, 0, 45); + auto b = + Manifold::Cube({1, 1, 1}, true).Rotate(0, 45, 0).Translate({2, 0, 0}); + + float distance = a.MinGap(b, 0.7f); + + EXPECT_FLOAT_EQ(distance, 2 - sqrt(2)); +} + +TEST(Manifold, MinGapClosestPointOnTriangleFace) { + auto a = Manifold::Cube(); + auto b = Manifold::Cube().Scale({10, 10, 10}).Translate({2, -5, -1}); + + float distance = a.MinGap(b, 1.1f); + + EXPECT_FLOAT_EQ(distance, 1); +} + +TEST(Manifold, MingapAfterTransformations) { + auto a = Manifold::Sphere(1, 512).Rotate(30, 30, 30); + auto b = + Manifold::Sphere(1, 512).Scale({3, 1, 1}).Rotate(0, 90, 45).Translate( + {3, 0, 0}); + + float distance = a.MinGap(b, 1.1f); + + ASSERT_NEAR(distance, 1, 0.001f); +} + +TEST(Manifold, MingapStretchyBracelet) { + auto a = StretchyBracelet(); + auto b = StretchyBracelet().Translate({0, 0, 20}); + + float distance = a.MinGap(b, 10); + + ASSERT_NEAR(distance, 5, 0.001f); +} + +TEST(Manifold, MinGapAfterTransformationsOutOfBounds) { + auto a = Manifold::Sphere(1, 512).Rotate(30, 30, 30); + auto b = + Manifold::Sphere(1, 512).Scale({3, 1, 1}).Rotate(0, 90, 45).Translate( + {3, 0, 0}); + + float distance = a.MinGap(b, 0.95f); + + ASSERT_NEAR(distance, 0.95f, 0.001f); +} +TEST(Manifold, TriangleDistanceClosestPointsOnVertices) { + std::array p = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0, 0}, + glm::vec3{0, 1, 0}}; + + std::array q = {glm::vec3{2, 0, 0}, glm::vec3{4, 0, 0}, + glm::vec3{3, 1, 0}}; + + float distance = DistanceTriangleTriangleSquared(p, q); + + EXPECT_FLOAT_EQ(distance, 1); +} + +TEST(Manifold, TriangleDistanceClosestPointOnEdge) { + std::array p = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0, 0}, + glm::vec3{0, 1, 0}}; + + std::array q = {glm::vec3{-1, 2, 0}, glm::vec3{1, 2, 0}, + glm::vec3{0, 3, 0}}; + + float distance = DistanceTriangleTriangleSquared(p, q); + + EXPECT_FLOAT_EQ(distance, 1); +} + +TEST(Manifold, TriangleDistanceClosestPointOnEdge2) { + std::array p = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0, 0}, + glm::vec3{0, 1, 0}}; + + std::array q = {glm::vec3{1, 1, 0}, glm::vec3{3, 1, 0}, + glm::vec3{2, 2, 0}}; + + float distance = DistanceTriangleTriangleSquared(p, q); + + EXPECT_FLOAT_EQ(distance, 0.5f); +} + +TEST(Manifold, TriangleDistanceClosestPointOnFace) { + std::array p = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0, 0}, + glm::vec3{0, 1, 0}}; + + std::array q = {glm::vec3{-1, 2, -0.5f}, glm::vec3{1, 2, -0.5f}, + glm::vec3{0, 2, 1.5f}}; + + float distance = DistanceTriangleTriangleSquared(p, q); + + EXPECT_FLOAT_EQ(distance, 1); +} + +TEST(Manifold, TriangleDistanceOverlapping) { + std::array p = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0, 0}, + glm::vec3{0, 1, 0}}; + + std::array q = {glm::vec3{-1, 0, 0}, glm::vec3{1, 0.5f, 0}, + glm::vec3{0, 1, 0}}; + + float distance = DistanceTriangleTriangleSquared(p, q); + + EXPECT_FLOAT_EQ(distance, 0); +}