Add Convex Decomposition and Minkowski Sum

.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

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
 

bindings/c/manifoldc.cpp

@@ -405,6 +405,25 @@ ManifoldManifoldVec *manifold_decompose(void *mem, ManifoldManifold *m) {
   return to_c(new (mem) std::vector<Manifold>(comps));
 }
 
+ManifoldManifold *manifold_fracture(void *mem, ManifoldManifold *m,
+                                    ManifoldVec3 *points, float *weights,
+                                    size_t length) {
+  std::vector<glm::dvec3> pts(length);
+  std::vector<double> 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<Manifold>(comps));
+}
+
+ManifoldManifold *manifold_convex_decomposition(void *mem,
+                                                ManifoldManifold *m) {
+  auto comps = from_c(m)->ConvexDecomposition();
+  return to_c(new (mem) std::vector<Manifold>(comps));
+}
+
 ManifoldMeshGL *manifold_get_meshgl(void *mem, ManifoldManifold *m) {
   auto mesh = from_c(m)->GetMeshGL();
   return to_c(new (mem) MeshGL(mesh));

bindings/python/manifold3d.cpp

@@ -289,6 +289,12 @@ NB_MODULE(manifold3d, m) {
           [](std::vector<glm::vec3> 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<nb::numpy, const double, nb::c_contig,
+                               nb::shape<nb::any, 3>> &pts,
+             const nb::ndarray<nb::numpy, const double, nb::c_contig,
+                               nb::shape<nb::any>> &weights) {
+            std::vector<glm::dvec3> pts_vec;
+            std::vector<double> 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"

bindings/wasm/bindings.cpp

@@ -152,6 +152,8 @@ EMSCRIPTEN_BINDINGS(whatever) {
       .function("_Mirror", &Manifold::Mirror)
       .function("_Decompose", select_overload<std::vector<Manifold>() const>(
                                   &Manifold::Decompose))
+      .function("_Fracture", &Manifold::Fracture)
+      .function("_ConvexDecomposition", &Manifold::ConvexDecomposition)
       .function("isEmpty", &Manifold::IsEmpty)
       .function("status", &Manifold::Status)
       .function("numVert", &Manifold::NumVert)

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 {

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()

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 <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/convex_decomposition_3.h>
+#include <CGAL/minkowski_sum_3.h>
+
+#include <chrono>
+#include <iostream>
+
+#include "manifold.h"
+
+using namespace manifold;
+
+using Kernel = CGAL::Epeck;
+using Polyhedron = CGAL::Polyhedron_3<Kernel>;
+using HalfedgeDS = Polyhedron::HalfedgeDS;
+using NefPolyhedron = CGAL::Nef_polyhedron_3<Kernel>;
+
+template <class HDS>
+class BuildFromManifold : public CGAL::Modifier_base<HDS> {
+ 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<HDS> 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<HalfedgeDS> 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<HalfedgeDS>(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;
+}

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})

src/manifold/include/manifold.h

@@ -249,6 +249,26 @@ class Manifold {
   Manifold Hull() const;
   static Manifold Hull(const std::vector<Manifold>& manifolds);
   static Manifold Hull(const std::vector<glm::vec3>& pts);
+  static std::vector<Manifold> BatchHull(
+      const std::vector<std::vector<Manifold>>& manifolds);
+  ///@}
+
+  /** @name Voronoi Functions
+   */
+  ///@{
+  std::vector<Manifold> Fracture(const std::vector<glm::dvec3>& pts,
+                                 const std::vector<double>& weights) const;
+  std::vector<Manifold> Fracture(const std::vector<glm::vec3>& pts,
+                                 const std::vector<float>& weights) const;
+  std::vector<int> Manifold::ReflexFaces(double tolerance = 1e-8) const;
+  std::vector<Manifold> ConvexDecomposition() const;
+  ///@}
+
+  /** @name Minkowski Functions
+   */
+  ///@{
+  static Manifold Minkowski(const Manifold& a, const Manifold& b,
+                            bool useNaive = false);
   ///@}
 
   /** @name Testing hooks

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<glm::dvec3> 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

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<glm::dvec3> verts, int face) const;
 
   // edge_op.cu
   void SimplifyTopology();