diff --git a/.gitmodules b/.gitmodules index 338f77fd8..b58e714f5 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "bindings/python/third_party/nanobind"] path = bindings/python/third_party/nanobind url = https://github.com/wjakob/nanobind +[submodule "src/third_party/voro"] + path = src/third_party/voro + url = https://github.com/chr1shr/voro diff --git a/bindings/c/include/manifoldc.h b/bindings/c/include/manifoldc.h index 8203405ab..221a09b25 100644 --- a/bindings/c/include/manifoldc.h +++ b/bindings/c/include/manifoldc.h @@ -166,6 +166,11 @@ ManifoldManifold *manifold_revolve(void *mem, ManifoldCrossSection *cs, ManifoldManifold *manifold_compose(void *mem, ManifoldManifoldVec *ms); ManifoldManifoldVec *manifold_decompose(void *mem, ManifoldManifold *m); ManifoldManifold *manifold_as_original(void *mem, ManifoldManifold *m); +ManifoldManifoldVec *manifold_fracture(void *mem, ManifoldManifold *m, + ManifoldVec3 *points, float *weights, + size_t length); +ManifoldManifoldVec *manifold_convex_decomposition(void *mem, + ManifoldManifold *m); // Manifold Info diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp index 5e1ea9652..67519432f 100644 --- a/bindings/c/manifoldc.cpp +++ b/bindings/c/manifoldc.cpp @@ -405,6 +405,25 @@ ManifoldManifoldVec *manifold_decompose(void *mem, ManifoldManifold *m) { return to_c(new (mem) std::vector(comps)); } +ManifoldManifold *manifold_fracture(void *mem, ManifoldManifold *m, + ManifoldVec3 *points, float *weights, + size_t length) { + std::vector pts(length); + std::vector wts(length); + for (int i = 0; i < length; ++i) { + pts[i] = {points[i].x, points[i].y, points[i].z}; + wts[i] = weights[i]; + } + auto comps = from_c(m)->Fracture(pts, wts); + return to_c(new (mem) std::vector(comps)); +} + +ManifoldManifold *manifold_convex_decomposition(void *mem, + ManifoldManifold *m) { + auto comps = from_c(m)->ConvexDecomposition(); + return to_c(new (mem) std::vector(comps)); +} + ManifoldMeshGL *manifold_get_meshgl(void *mem, ManifoldManifold *m) { auto mesh = from_c(m)->GetMeshGL(); return to_c(new (mem) MeshGL(mesh)); diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp index 6e0bd1357..27f04846f 100644 --- a/bindings/python/manifold3d.cpp +++ b/bindings/python/manifold3d.cpp @@ -289,6 +289,12 @@ NB_MODULE(manifold3d, m) { [](std::vector pts) { return Manifold::Hull(pts); }, nb::arg("pts"), "Compute the convex hull enveloping a set of 3d points.") + .def_static( + "batch_batch_hull", Manifold::BatchHull, + "Compute the convex hulls enveloping multiple sets of manifolds.") + .def_static("minkowski", Manifold::Minkowski, nb::arg("a"), nb::arg("b"), + nb::arg("useNaive"), + "Compute the minkowski sum of two manifolds.") .def( "transform", &Manifold::Transform, nb::arg("m"), "Transform this Manifold in space. The first three columns form a " @@ -528,6 +534,32 @@ NB_MODULE(manifold3d, m) { .def("project", &Manifold::Project, "Returns a cross section representing the projected outline of this " "object onto the X-Y plane.") + .def( + "fracture", + [](Manifold &self, + const nb::ndarray> &pts, + const nb::ndarray> &weights) { + std::vector pts_vec; + std::vector weights_vec; + auto pointData = pts.data(); + auto weightData = weights.data(); + for (int i = 0; i < pts.shape(0) * pts.shape(1); + i += pts.shape(1)) { + pts_vec.push_back( + {pointData[i], pointData[i + 1], pointData[i + 2]}); + } + for (int i = 0; i < weights.shape(0); i++) { + weights_vec.push_back(weightData[i]); + } + return self.Fracture(pts_vec, weights_vec); + }, + "This operation computes the fracturing of this Manifold into the " + "chunks around the supplied points.") + .def("convex_decomposition", &Manifold::ConvexDecomposition, + "This operation computes the fracturing of this Manifold into the " + "minimal number of representative convex pieces.") .def("status", &Manifold::Status, "Returns the reason for an input Mesh producing an empty Manifold. " "This Status only applies to Manifolds newly-created from an input " @@ -622,6 +654,11 @@ NB_MODULE(manifold3d, m) { ":param center: Set to true to shift the center to the origin. " "Default is origin at the bottom.") .def_static( + "sphere", + [](float radius, + int circularSegments = + 0) { return Manifold::Sphere(radius, circularSegments); }, + nb::arg("radius"), nb::arg("circular_segments") = 0, "sphere", &Manifold::Sphere, nb::arg("radius"), nb::arg("circular_segments") = 0, "Constructs a geodesic sphere of a given radius.\n" diff --git a/bindings/wasm/bindings.cpp b/bindings/wasm/bindings.cpp index 778b94938..b7ada08aa 100644 --- a/bindings/wasm/bindings.cpp +++ b/bindings/wasm/bindings.cpp @@ -152,6 +152,8 @@ EMSCRIPTEN_BINDINGS(whatever) { .function("_Mirror", &Manifold::Mirror) .function("_Decompose", select_overload() const>( &Manifold::Decompose)) + .function("_Fracture", &Manifold::Fracture) + .function("_ConvexDecomposition", &Manifold::ConvexDecomposition) .function("isEmpty", &Manifold::IsEmpty) .function("status", &Manifold::Status) .function("numVert", &Manifold::NumVert) diff --git a/bindings/wasm/bindings.js b/bindings/wasm/bindings.js index ac49132c5..b481128af 100644 --- a/bindings/wasm/bindings.js +++ b/bindings/wasm/bindings.js @@ -289,6 +289,52 @@ Module.setup = function() { return result; }; + Module.Manifold.prototype.fracture = function(points, weights) { + if (points instanceof Array) { + let pts = new Module.Vector_vec3(); + let wts = new Module.Vector_f32(); + for (const m of points) { + if (m instanceof Array && typeof m[0] == 'number' && m.length >= 3) { + pts.push_back({x: m[0], y: m[1], z: m[2]}); + if (m.length >= 4) { + wts.push_back(m[3]); + } + } else if (m.x) { + pts.push_back(m); + } else { + pushVec3(pts, m); + } + } + if (weights instanceof Array && weights.length == points.length) { + for (const m of weights) { + if (typeof m == 'number') { + wts.push_back(m); + } + } + } else { + if (wts.size() == 0) { + for (const m of points) { + wts.push_back(1.0); + } + } + } + + const vec = this._Fracture(pts, wts); + pts.delete(); + wts.delete(); + const result = fromVec(vec); + vec.delete(); + return result; + } + }; + + Module.Manifold.prototype.convexDecomposition = function() { + const vec = this._ConvexDecomposition(); + const result = fromVec(vec); + vec.delete(); + return result; + }; + Module.Manifold.prototype.boundingBox = function() { const result = this._boundingBox(); return { diff --git a/extras/CMakeLists.txt b/extras/CMakeLists.txt index 27824f270..d23ca7e6a 100644 --- a/extras/CMakeLists.txt +++ b/extras/CMakeLists.txt @@ -39,11 +39,18 @@ if(MANIFOLD_EXPORT) endif() if(BUILD_TEST_CGAL) - add_executable(perfTestCGAL perf_test_cgal.cpp) find_package(CGAL REQUIRED COMPONENTS Core) find_package(Boost REQUIRED COMPONENTS thread) + + add_executable(perfTestCGAL perf_test_cgal.cpp) target_compile_definitions(perfTestCGAL PRIVATE CGAL_USE_GMPXX) target_link_libraries(perfTestCGAL manifold CGAL::CGAL CGAL::CGAL_Core Boost::thread) target_compile_options(perfTestCGAL PRIVATE ${MANIFOLD_FLAGS}) target_compile_features(perfTestCGAL PUBLIC cxx_std_17) + + add_executable(perfTestCGALCD perf_test_cgal_cd.cpp) + target_compile_definitions(perfTestCGALCD PRIVATE CGAL_USE_GMPXX) + target_link_libraries(perfTestCGALCD manifold CGAL::CGAL CGAL::CGAL_Core Boost::thread) + target_compile_options(perfTestCGALCD PRIVATE ${MANIFOLD_FLAGS}) + target_compile_features(perfTestCGALCD PUBLIC cxx_std_17) endif() diff --git a/extras/perf_test_cgal_cd.cpp b/extras/perf_test_cgal_cd.cpp new file mode 100644 index 000000000..3f7eb32d3 --- /dev/null +++ b/extras/perf_test_cgal_cd.cpp @@ -0,0 +1,120 @@ +// Copyright 2023 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. + +#include +#include +#include +#include +#include + +#include +#include + +#include "manifold.h" + +using namespace manifold; + +using Kernel = CGAL::Epeck; +using Polyhedron = CGAL::Polyhedron_3; +using HalfedgeDS = Polyhedron::HalfedgeDS; +using NefPolyhedron = CGAL::Nef_polyhedron_3; + +template +class BuildFromManifold : public CGAL::Modifier_base { + public: + using Vertex = typename HDS::Vertex; + using Point = typename Vertex::Point; + BuildFromManifold(const Manifold manifold) : manifold(manifold) {} + void operator()(HDS &hds) { + // Postcondition: hds is a valid polyhedral surface. + CGAL::Polyhedron_incremental_builder_3 B(hds, true); + auto mesh = manifold.GetMeshGL(); + B.begin_surface(mesh.NumVert(), mesh.NumTri(), mesh.NumTri() * 3); + + for (size_t i = 0; i < mesh.vertProperties.size(); i += mesh.numProp) { + B.add_vertex(Point(mesh.vertProperties[i], mesh.vertProperties[i + 1], + mesh.vertProperties[i + 2])); + } + for (size_t i = 0; i < mesh.triVerts.size(); i += 3) { + B.begin_facet(); + for (const int j : {0, 1, 2}) B.add_vertex_to_facet(mesh.triVerts[i + j]); + B.end_facet(); + } + B.end_surface(); + } + + private: + const Manifold manifold; +}; + +int main(int argc, char **argv) { + Manifold spherecube = + Manifold::Cube(glm::vec3(1), true) - Manifold::Sphere(0.6, 100); + Manifold smallsphere = Manifold::Sphere(0.1, 20); + + BuildFromManifold build(spherecube); + std::cout << "nTri = " << spherecube.NumTri() << std::endl; + + auto start = std::chrono::high_resolution_clock::now(); + Polyhedron poly; + poly.delegate(build); + std::cout << "to Polyhedron took " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + start = std::chrono::high_resolution_clock::now(); + NefPolyhedron np(poly); + std::cout << "conversion to Nef Polyhedron took " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + start = std::chrono::high_resolution_clock::now(); + auto convexDecomposition = spherecube.ConvexDecomposition(); + std::cout << "[MANIFOLD] decomposed into " << convexDecomposition.size() + << " parts in " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + start = std::chrono::high_resolution_clock::now(); + auto generalMinkowskiSum = Manifold::Minkowski(spherecube, smallsphere); + std::cout << "[MANIFOLD] general minkowski summed in " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + start = std::chrono::high_resolution_clock::now(); + auto naiveMinkowskiSum = Manifold::Minkowski(spherecube, smallsphere, true); + std::cout << "[MANIFOLD] naive minkowski summed in " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + start = std::chrono::high_resolution_clock::now(); + CGAL::convex_decomposition_3(np); + std::cout << "[CGAL] decomposed into " << np.number_of_volumes() + << " parts in " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + // Create the Small Sphere NEF Polyhedron for Minkowski Summing + Polyhedron poly2; + poly.delegate(BuildFromManifold(smallsphere)); + NefPolyhedron np2(poly); + + start = std::chrono::high_resolution_clock::now(); + CGAL::minkowski_sum_3(np, np2); + std::cout << "[CGAL] minkowski summed in " + << (std::chrono::high_resolution_clock::now() - start).count() / 1e9 + << " sec" << std::endl; + + return 0; +} diff --git a/src/manifold/CMakeLists.txt b/src/manifold/CMakeLists.txt index 72cff48d7..e7a5749b7 100644 --- a/src/manifold/CMakeLists.txt +++ b/src/manifold/CMakeLists.txt @@ -20,7 +20,7 @@ add_library(${PROJECT_NAME} ${SOURCE_FILES}) target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR}/include) target_link_libraries(${PROJECT_NAME} PUBLIC utilities cross_section - PRIVATE collider polygon ${MANIFOLD_INCLUDE} Clipper2 quickhull + PRIVATE collider polygon ${MANIFOLD_INCLUDE} Clipper2 quickhull voro++ ) target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS}) diff --git a/src/manifold/include/manifold.h b/src/manifold/include/manifold.h index 11366d5c8..b7798f5c6 100644 --- a/src/manifold/include/manifold.h +++ b/src/manifold/include/manifold.h @@ -249,6 +249,26 @@ class Manifold { Manifold Hull() const; static Manifold Hull(const std::vector& manifolds); static Manifold Hull(const std::vector& pts); + static std::vector BatchHull( + const std::vector>& manifolds); + ///@} + + /** @name Voronoi Functions + */ + ///@{ + std::vector Fracture(const std::vector& pts, + const std::vector& weights) const; + std::vector Fracture(const std::vector& pts, + const std::vector& weights) const; + std::vector Manifold::ReflexFaces(double tolerance = 1e-8) const; + std::vector ConvexDecomposition() const; + ///@} + + /** @name Minkowski Functions + */ + ///@{ + static Manifold Minkowski(const Manifold& a, const Manifold& b, + bool useNaive = false); ///@} /** @name Testing hooks diff --git a/src/manifold/src/face_op.cpp b/src/manifold/src/face_op.cpp index 331287655..3ba371868 100644 --- a/src/manifold/src/face_op.cpp +++ b/src/manifold/src/face_op.cpp @@ -316,4 +316,33 @@ CrossSection Manifold::Impl::Project() const { return CrossSection(polys).Simplify(precision_); } + +glm::dvec4 Manifold::Impl::Circumcircle(Vec verts, int face) const { + glm::dvec3 va = verts[this->halfedge_[(face * 3) + 0].startVert]; + glm::dvec3 vb = verts[this->halfedge_[(face * 3) + 1].startVert]; + glm::dvec3 vc = verts[this->halfedge_[(face * 3) + 2].startVert]; + + glm::dvec3 a = va - vc; + glm::dvec3 b = vb - vc; + glm::dvec3 c = va - vb; + double a_length = glm::length(a); + double b_length = glm::length(b); + double c_length = glm::length(c); + glm::dvec3 numerator = + glm::cross((((a_length * a_length) * b) - ((b_length * b_length) * a)), + glm::cross(a, b)); + double crs = glm::length(glm::cross(a, b)); + double denominator = 2.0 * (crs * crs); + glm::dvec3 circumcenter = (numerator / denominator) + vc; + double circumradius = glm::length(circumcenter - vc); + + double max_length = std::fmax(a_length, std::fmax(b_length, c_length)); + double min_length = std::fmin(a_length, std::fmin(b_length, c_length)); + if (max_length / min_length > 15.0) { + circumradius *= -1.0; // Mark this triangle as degenerate + } + + return glm::dvec4(circumcenter.x, circumcenter.y, circumcenter.z, + circumradius); +} } // namespace manifold diff --git a/src/manifold/src/impl.h b/src/manifold/src/impl.h index 36e37fcc0..3889d4dcb 100644 --- a/src/manifold/src/impl.h +++ b/src/manifold/src/impl.h @@ -121,6 +121,7 @@ struct Manifold::Impl { glm::mat3x2 projection) const; CrossSection Slice(float height) const; CrossSection Project() const; + glm::dvec4 Circumcircle(Vec verts, int face) const; // edge_op.cu void SimplifyTopology(); diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 841339aac..1d612ed04 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -15,12 +15,15 @@ #include #include #include +#include +#include #include "QuickHull.hpp" #include "boolean3.h" #include "csg_tree.h" #include "impl.h" #include "par.h" +#include "voro++.hh" namespace { using namespace manifold; @@ -61,6 +64,68 @@ struct UpdateProperties { } }; +struct ComputeVoronoiCell { + voro::container_poly* container; + const Manifold* original; + std::vector* output; + void operator()(thrust::tuple inOut) { + const int idx = thrust::get<0>(inOut); + glm::ivec3& cell_idx = thrust::get<1>(inOut); + glm::dvec4& cell_pos = thrust::get<2>(inOut); + voro::voronoicell_neighbor c(*container); + + if (container->compute_cell(c, cell_idx.y, cell_idx.z)) { + std::vector verts; + verts.reserve(c.p); + for (size_t i = 0; i < c.p; i++) { + verts.push_back(glm::vec3(cell_pos.x + 0.5 * c.pts[(4 * i) + 0], + cell_pos.y + 0.5 * c.pts[(4 * i) + 1], + cell_pos.z + 0.5 * c.pts[(4 * i) + 2])); + } + (*output)[cell_idx.x] = Manifold::Hull(verts) ^ *original; + verts.clear(); + } else { + std::cout << "[ERROR] Degenerate Voronoi Cell at Index: " << cell_idx.x + << std::endl; + } + } +}; + +struct ComputeTriangleHull { + const Manifold* bM; + std::vector* vertPos; + std::vector* output; + void operator()(thrust::tuple inOut) { + const int idx = thrust::get<0>(inOut); + glm::ivec3& tri = thrust::get<1>(inOut); + (*output)[idx] = Manifold::Hull({(*bM).Translate((*vertPos)[tri.x]), + (*bM).Translate((*vertPos)[tri.y]), + (*bM).Translate((*vertPos)[tri.z])}); + } +}; + +struct ComputeTriangleTriangleHull { + std::vector* vertPosAPtr; + std::vector* vertPosBPtr; + std::vector* output; + void operator()(thrust::tuple> inOut) { + const int idx = thrust::get<0>(inOut); + std::pair tris = thrust::get<1>(inOut); + auto vertPosA = *vertPosAPtr; + auto vertPosB = *vertPosBPtr; + (*output)[idx] = + Manifold::Hull({vertPosA[tris.first.x] + vertPosB[tris.second.x], + vertPosA[tris.first.x] + vertPosB[tris.second.y], + vertPosA[tris.first.x] + vertPosB[tris.second.z], + vertPosA[tris.first.y] + vertPosB[tris.second.x], + vertPosA[tris.first.y] + vertPosB[tris.second.y], + vertPosA[tris.first.y] + vertPosB[tris.second.z], + vertPosA[tris.first.z] + vertPosB[tris.second.x], + vertPosA[tris.first.z] + vertPosB[tris.second.y], + vertPosA[tris.first.z] + vertPosB[tris.second.z]}); + } +}; + Manifold Halfspace(Box bBox, glm::vec3 normal, float originOffset) { normal = glm::normalize(normal); Manifold cutter = @@ -859,4 +924,280 @@ Manifold Manifold::Hull() const { return Hull(GetMesh().vertPos); } Manifold Manifold::Hull(const std::vector& manifolds) { return Compose(manifolds).Hull(); } + +/** + * Compute the convex hulls enveloping many sets of manifolds. + * + * @param manifolds A vector of vectors of manifolds over which compute hulls. + */ +std::vector Manifold::BatchHull( + const std::vector>& manifolds) { + ZoneScoped; + std::vector output; + output.reserve(manifolds.size()); + output.resize(manifolds.size()); + thrust::for_each_n( + thrust::host, zip(manifolds.begin(), output.begin()), manifolds.size(), + [](thrust::tuple, Manifold&> inOut) { + thrust::get<1>(inOut) = Manifold::Hull(thrust::get<0>(inOut)); + }); + return output; +} + +/** + * Compute the voronoi fracturing of this Manifold into convex chunks. + * + * @param pts A vector of points over which to fracture the manifold. + * @param wts A vector of weights controlling the relative size of each chunk. + */ +std::vector Manifold::Fracture(const std::vector& pts, + const std::vector& wts) const { + ZoneScoped; + std::vector output; + output.reserve(pts.size()); + output.resize(pts.size()); + + Box bounds = BoundingBox(); + glm::vec3 min = bounds.min - 0.1f; + glm::vec3 max = bounds.max + 0.1f; + double V = (max.x - min.x) * (max.y - min.y) * (max.z - min.z); + double Nthird = powf((double)pts.size() / V, 1.0f / 3.0f); + voro::container_poly container(min.x, max.x, min.y, max.y, min.z, max.z, + std::round(Nthird * (max.x - min.x)), + std::round(Nthird * (max.y - min.y)), + std::round(Nthird * (max.z - min.z)), // + false, false, false, pts.size()); + + bool hasWeights = wts.size() == pts.size(); + for (size_t i = 0; i < pts.size(); i++) { + container.put(i, pts[i].x, pts[i].y, pts[i].z, hasWeights ? wts[i] : 1.0f); + } + + // Prepare Parallel Voronoi Computation + std::vector cellIndices; + std::vector cellPosWeight; + voro::c_loop_all vl(container); + if (vl.start()) do { + int id; + double x, y, z, r; + vl.pos(id, x, y, z, r); + cellIndices.push_back({id, vl.ijk, vl.q}); + cellPosWeight.push_back({x, y, z, r}); + } while (vl.inc()); + + ComputeVoronoiCell computeVoronoiCell; + computeVoronoiCell.container = &container; + computeVoronoiCell.original = this; + computeVoronoiCell.output = &output; + thrust::for_each_n( + thrust::host, zip(countAt(0), cellIndices.begin(), cellPosWeight.begin()), + cellIndices.size(), computeVoronoiCell); + return output; +} +std::vector Manifold::Fracture( + const std::vector& pts, + const std::vector& weights) const { + std::vector highpVerts(pts.size()); + std::vector highpWeights(weights.size()); + for (size_t i = 0; i < pts.size(); i++) { + highpVerts[i] = pts[i]; + highpWeights[i] = weights[i]; + } + return Fracture(highpVerts, highpWeights); +} + +std::vector Manifold::ReflexFaces(double tolerance) const { + std::unordered_set uniqueReflexFaceSet; + const Impl& impl = *GetCsgLeafNode().GetImpl(); + for (size_t i = 0; i < impl.halfedge_.size(); i++) { + Halfedge halfedge = impl.halfedge_[i]; + int faceA = halfedge.face; + int faceB = impl.halfedge_[halfedge.pairedHalfedge].face; + glm::dvec3 tangent = + glm::cross((glm::dvec3)impl.faceNormal_[faceA], + (glm::dvec3)impl.vertPos_[impl.halfedge_[i].endVert] - + (glm::dvec3)impl.vertPos_[impl.halfedge_[i].startVert]); + double tangentProjection = + glm::dot((glm::dvec3)impl.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 uniqueFaces; // Copy to a vector for indexed access + uniqueFaces.insert(uniqueFaces.end(), uniqueReflexFaceSet.begin(), + uniqueReflexFaceSet.end()); + return uniqueFaces; +} + +std::vector Manifold::ConvexDecomposition() const { + ZoneScoped; + + // Early-Exit if the manifold is already convex + std::vector uniqueFaces = ReflexFaces(); + if (uniqueFaces.size() == 0) { + return std::vector(1, *this); + } + + const Impl& impl = *GetCsgLeafNode().GetImpl(); + std::vector circumcenters(uniqueFaces.size()); + std::vector circumradii(uniqueFaces.size()); + Vec dVerts(impl.vertPos_.size()); + for (size_t i = 0; i < impl.vertPos_.size(); i++) { + dVerts[i] = impl.vertPos_[i]; + } + std::vector debugShapes; + for (size_t i = 0; i < uniqueFaces.size(); i++) { + glm::dvec4 circumcircle = impl.Circumcircle(dVerts, uniqueFaces[i]); + circumcenters[i] = + glm::dvec3(circumcircle.x, circumcircle.y, circumcircle.z); + circumradii[i] = circumcircle.w; + + // Debug Draw the Circumcenter and Triangle of Degenerate Triangles + // if (circumcircle.w < 0.0) { + // std::cout << "Circumradius: " << circumcircle.w << std::endl; + // debugShapes.push_back(Hull( + // {Manifold::Sphere(circumcircle.w * 0.1, 10) + // .Translate(glm::vec3( + // dVerts[impl.halfedge_[(uniqueFaces[i] * 3) + + // 0].startVert])), + // Manifold::Sphere(circumcircle.w * 0.1, 10) + // .Translate(glm::vec3( + // dVerts[impl.halfedge_[(uniqueFaces[i] * 3) + + // 1].startVert])), + // Manifold::Sphere(circumcircle.w * 0.1, 10) + // .Translate(glm::vec3( + // dVerts[impl.halfedge_[(uniqueFaces[i] * 3) + + // 2].startVert])) + // })); + // + // debugShapes.push_back(Manifold::Sphere(circumcircle.w * 0.2, 10) + // .Translate(glm::vec3(circumcenters[i]))); + // } + } + + for (size_t i = 0; i < circumcenters.size() - 1; i++) { + for (size_t j = circumcenters.size() - 1; j > i; j--) { + if (glm::distance(circumcenters[i], circumcenters[j]) < 1e-9 || + circumradii[j] < 0.0) { + std::vector::iterator it = circumcenters.begin(); + std::advance(it, j); + circumcenters.erase(it); + std::vector::iterator it2 = circumradii.begin(); + std::advance(it2, j); + circumradii.erase(it2); + } + } + } + + std::vector output = Fracture(circumcenters, circumradii); + output.insert(output.end(), debugShapes.begin(), debugShapes.end()); + return output; +} + +/** + * Compute the minkowski sum of two manifolds. + * + * @param a The first manifold in the sum. + * @param b The second manifold in the sum. + */ +Manifold Manifold::Minkowski(const Manifold& a, const Manifold& b, + bool useNaive) { + std::vector> composedParts; + std::vector composedHulls({a}); + if (!useNaive) { // Use the general method + std::vector aDecomposition = a.ConvexDecomposition(); + std::vector bDecomposition = b.ConvexDecomposition(); + for (Manifold aPart : aDecomposition) { + manifold::Mesh aMesh = aPart.GetMesh(); + for (Manifold bPart : bDecomposition) { + std::vector abComposition; + for (glm::vec3 vertexPosition : aMesh.vertPos) { + abComposition.push_back(bPart.Translate(vertexPosition)); + } + composedParts.push_back(abComposition); + } + } + } else { // Use the naive method + bool aConvex = a.ReflexFaces().size() == 0; + bool bConvex = b.ReflexFaces().size() == 0; + + // If the convex manifold was supplied first, swap them! + Manifold aM = a, bM = b; + if (aConvex && !bConvex) { + aM = b; + bM = a; + aConvex = !aConvex; + bConvex = !bConvex; + } + + manifold::Mesh aMesh = aM.GetMesh(); + + // Convex-Convex Minkowski: Very Fast + if (aConvex && bConvex) { + std::vector simpleHull; + for (glm::vec3 vertex : aMesh.vertPos) { + simpleHull.push_back({bM.Translate(vertex)}); + } + composedHulls.push_back(Manifold::Hull(simpleHull)); + // Convex - Non-Convex Minkowski: Slower + } else if (!aConvex && bConvex) { + // Speed Equivalent? + // for (glm::ivec3 vertexIndices : aMesh.triVerts) { + // composedParts.push_back({bM.Translate(aMesh.vertPos[vertexIndices.x]), + // bM.Translate(aMesh.vertPos[vertexIndices.y]), + // bM.Translate(aMesh.vertPos[vertexIndices.z])}); + //} + composedHulls.resize(aMesh.triVerts.size() + 1); + thrust::for_each_n( + thrust::host, zip(countAt(1), aMesh.triVerts.begin()), + aMesh.triVerts.size(), + ComputeTriangleHull({&bM, &aMesh.vertPos, &composedHulls})); + // Non-Convex - Non-Convex Minkowski: Very Slow + } else if (!aConvex && !bConvex) { + manifold::Mesh bMesh = bM.GetMesh(); + + // Speed Equivalent? + // for (glm::ivec3 aVertexIndices : aMesh.triVerts) { + // for (glm::ivec3 bVertexIndices : bMesh.triVerts) { + // composedHulls.push_back(Manifold::Hull( + // {aMesh.vertPos[aVertexIndices.x] + + // bMesh.vertPos[bVertexIndices.x], + // aMesh.vertPos[aVertexIndices.x] + + // bMesh.vertPos[bVertexIndices.y], + // aMesh.vertPos[aVertexIndices.x] + + // bMesh.vertPos[bVertexIndices.z], + // aMesh.vertPos[aVertexIndices.y] + + // bMesh.vertPos[bVertexIndices.x], + // aMesh.vertPos[aVertexIndices.y] + + // bMesh.vertPos[bVertexIndices.y], + // aMesh.vertPos[aVertexIndices.y] + + // bMesh.vertPos[bVertexIndices.z], + // aMesh.vertPos[aVertexIndices.z] + + // bMesh.vertPos[bVertexIndices.x], + // aMesh.vertPos[aVertexIndices.z] + + // bMesh.vertPos[bVertexIndices.y], + // aMesh.vertPos[aVertexIndices.z] + + // bMesh.vertPos[bVertexIndices.z]})); + // } + //} + + std::vector> trianglePairs; + for (glm::ivec3 aVertexIndices : aMesh.triVerts) { + for (glm::ivec3 bVertexIndices : bMesh.triVerts) { + trianglePairs.push_back({aVertexIndices, bVertexIndices}); + } + } + composedHulls.resize(trianglePairs.size() + 1); + thrust::for_each_n(thrust::host, zip(countAt(1), trianglePairs.begin()), + trianglePairs.size(), + ComputeTriangleTriangleHull( + {&aMesh.vertPos, &bMesh.vertPos, &composedHulls})); + } + } + auto newHulls = Manifold::BatchHull(composedParts); + composedHulls.insert(composedHulls.end(), newHulls.begin(), newHulls.end()); + return Manifold::BatchBoolean(composedHulls, manifold::OpType::Add); +} } // namespace manifold diff --git a/src/third_party/CMakeLists.txt b/src/third_party/CMakeLists.txt index 6ef64e832..91f17c4f0 100644 --- a/src/third_party/CMakeLists.txt +++ b/src/third_party/CMakeLists.txt @@ -7,3 +7,9 @@ add_subdirectory(clipper2/CPP) add_library(quickhull quickhull/QuickHull.cpp) target_include_directories(quickhull PUBLIC quickhull) target_compile_features(quickhull PUBLIC cxx_std_17) + +set(VORO_BUILD_SHARED_LIBS OFF CACHE BOOL "") +set(VORO_BUILD_EXAMPLES OFF CACHE BOOL "") +set(VORO_BUILD_CMD_LINE OFF CACHE BOOL "") +set(VORO_ENABLE_DOXYGEN OFF CACHE BOOL "") +add_subdirectory(voro) diff --git a/src/third_party/voro b/src/third_party/voro new file mode 160000 index 000000000..56d619faf --- /dev/null +++ b/src/third_party/voro @@ -0,0 +1 @@ +Subproject commit 56d619faf3479313399516ad71c32773c29be859 diff --git a/test/manifold_test.cpp b/test/manifold_test.cpp index 14884af8a..ed253e271 100644 --- a/test/manifold_test.cpp +++ b/test/manifold_test.cpp @@ -784,3 +784,21 @@ TEST(Manifold, EmptyHull) { {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0}}; EXPECT_TRUE(Manifold::Hull(coplanar).IsEmpty()); } + +TEST(Manifold, ConvexDecomposition) { + Manifold sphere = Manifold::Sphere(0.6, 20); + Manifold cube = Manifold::Cube({1.0, 1.0, 1.0}, true); + Manifold nonConvex = cube - sphere; + std::vector convexParts = nonConvex.ConvexDecomposition(); + + EXPECT_EQ(convexParts.size(), 224); + + float originalVolume = nonConvex.GetProperties().volume; + float convex_volume = 0.0; + Manifold manifold_union = convexParts[0].AsOriginal(); + for (Manifold cur_manifold : convexParts) { + manifold_union += cur_manifold.Hull(); + } + float union_volume = manifold_union.GetProperties().volume; + EXPECT_NEAR(originalVolume, union_volume, 1e-6); +}