Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Convex Decomposition and Minkowski Sum #663

Closed
wants to merge 37 commits into from
Closed
Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
37 commits
Select commit Hold shift + click to select a range
ef30739
Add the Voro++ Library
zalo Dec 16, 2023
56e026b
Add Basic Python Binding
zalo Dec 16, 2023
804d29a
Get Fracture Operation Working
zalo Dec 16, 2023
15579ef
Add Broken VoACD Implementation
zalo Dec 16, 2023
c039b8a
Fix late night typos and formatting
zalo Dec 16, 2023
3f6248d
Get Basic VoCD Working
zalo Dec 17, 2023
a32be23
Formatting
zalo Dec 17, 2023
2aef6ff
Formatting the Sequel
zalo Dec 17, 2023
21f1daa
Formatting Part 3
zalo Dec 17, 2023
f0b9faf
Final Formatting
zalo Dec 17, 2023
ae5114d
After Final Formatting
zalo Dec 17, 2023
595c933
Use the common form of glm::dvec
zalo Dec 18, 2023
395e94b
Also dvec4
zalo Dec 18, 2023
7e1fa10
Reformat again
zalo Dec 18, 2023
dfe8e62
Remove the non-submodule folder
zalo Dec 18, 2023
1034970
Switch Voro to a Git Submodule
zalo Dec 18, 2023
76ca0d9
Add a test
zalo Dec 18, 2023
4921354
Add untested bindings
zalo Dec 18, 2023
d8fd65d
Minor Formatting
zalo Dec 18, 2023
04ae42e
More formatting
zalo Dec 18, 2023
7909f59
Dummy Commit
zalo Dec 18, 2023
2505b15
Parallelize Voronoi Computation
zalo Dec 18, 2023
f86adc5
Add Early Exit
zalo Dec 18, 2023
c574580
Add a function for computing many convex hulls in parallel
zalo Dec 18, 2023
1cb3a83
Add Minkowski Function
zalo Dec 18, 2023
e84e05c
Fix Formatting
zalo Dec 18, 2023
a4d20c0
Remove Extra Qualification
zalo Dec 18, 2023
9c04c6d
added cgal convex decomposition test
pca006132 Dec 19, 2023
33c69d6
Merge branch 'master' into feat-voacd
zalo Dec 19, 2023
b6978bf
Update CGAL Test with Direct Comparison...
zalo Dec 19, 2023
7ab872b
Add Minkowski Sum Test Too
zalo Dec 19, 2023
3446bb0
Remove the need for Joggling
zalo Dec 19, 2023
258612c
Add Non-Convex - Non-Convex Implementation
zalo Dec 19, 2023
3dcd3a5
Handle alternate order submission
zalo Dec 19, 2023
6a8a407
Finish Multithreading Naive Function
zalo Dec 19, 2023
c4537e7
Fix Formatting
zalo Dec 19, 2023
19c58ca
Detect Degenerate Triangles
zalo Dec 20, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitmodules
Original file line number Diff line number Diff line change
Expand Up @@ -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
26 changes: 26 additions & 0 deletions bindings/python/manifold3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -524,6 +524,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 "
Expand Down
2 changes: 1 addition & 1 deletion src/manifold/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand Down
10 changes: 10 additions & 0 deletions src/manifold/include/manifold.h
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,16 @@ class Manifold {
static Manifold Hull(const std::vector<glm::vec3>& pts);
///@}

/** @name Voronoi Functions
*/
///@{
std::vector<Manifold> Fracture(const std::vector<glm::dvec3>& pts,
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is fracture really meant to be a public function, or is it just a means to ConvexDecomposition?

Copy link
Contributor Author

@zalo zalo Dec 18, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's worth making it public because it's a non-trival operation that generalizes splitting a singular cutting plane by limiting the influence of each half-space. This can help make modelling cuts and edits more "local". Like, if you wanted to take a complicated, winding assembly, and break it apart into multiple pieces, it's nice to be able to make cuts with a localized area of effect to avoid mangling the model far from the cut.

Also, it's useful for people looking to implement their own approximate decompositions based on their particular speed and accuracy needs.

And it's good for modelling lattices and packings, which is one of the primary things Voro++ is used for... it could even help replicate NTopology lattice generation functionality if each cell is contracted and then subtracted from a solid region.

There's a universe with a dual version of this function, where the user directly specifies splitting planes and their relative area of influence. This could be implemented by inserting two voronoi cells right next to each other with the plane normal as the offset between them. Perhaps this would be more intuitive?

const std::vector<double>& weights) const;
std::vector<Manifold> Fracture(const std::vector<glm::vec3>& pts,
const std::vector<float>& weights) const;
std::vector<Manifold> ConvexDecomposition() const;
///@}

/** @name Testing hooks
* These are just for internal testing.
*/
Expand Down
20 changes: 20 additions & 0 deletions src/manifold/src/face_op.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -316,4 +316,24 @@ 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];
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about a mat3 and a for loop?


glm::dvec3 a = va - vc;
glm::dvec3 b = vb - vc;
double a_length = glm::length(a);
double b_length = glm::length(b);
glm::dvec3 numerator =
glm::cross((((a_length * a_length) * b) - ((b_length * b_length) * a)),
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we remove a few parens?

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);
return glm::dvec4(circumcenter.x, circumcenter.y, circumcenter.z,
circumradius);
}
} // namespace manifold
1 change: 1 addition & 0 deletions src/manifold/src/impl.h
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,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();
Expand Down
140 changes: 140 additions & 0 deletions src/manifold/src/manifold.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,15 @@
#include <algorithm>
#include <map>
#include <numeric>
#include <random>
#include <unordered_set>

#include "QuickHull.hpp"
#include "boolean3.h"
#include "csg_tree.h"
#include "impl.h"
#include "par.h"
#include "voro++.hh"

namespace {
using namespace manifold;
Expand Down Expand Up @@ -845,4 +848,141 @@ Manifold Manifold::Hull() const { return Hull(GetMesh().vertPos); }
Manifold Manifold::Hull(const std::vector<Manifold>& manifolds) {
return Compose(manifolds).Hull();
}

// TODO: Handle Joggling in the Fracture Function Directly?
/**
* Compute the voronoi fracturing of this Manifold into convex chunks.
*
* @param pts A vector of points over which to fracture the manifold.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we have some detail on what these inputs do and the purpose of this function?

* @param pts A vector of weights controlling the relative size of each chunk.
*/
std::vector<Manifold> Manifold::Fracture(
const std::vector<glm::dvec3>& pts,
const std::vector<double>& weights) const {
std::vector<Manifold> output;
output.reserve(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 = weights.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 ? weights[i] : 1.0f);
}
voro::voronoicell c(container);
voro::c_loop_all vl(container);
if (vl.start()) do { // TODO: Parallelize this loop!
if (container.compute_cell(c, vl)) {
std::vector<glm::vec3> verts;
verts.reserve(c.p);
int id;
double x, y, z, r;
vl.pos(id, x, y, z, r);
for (size_t i = 0; i < c.p; i++) {
verts.push_back(glm::vec3(x + 0.5 * c.pts[(vl.ps * i) + 0],
y + 0.5 * c.pts[(vl.ps * i) + 1],
z + 0.5 * c.pts[(vl.ps * i) + 2]));
}
Manifold outputManifold =
Hull(verts) ^
*this; // TODO: Replace this intersection with a voro++ wall
// implementation or cutting the cell directly...
if (outputManifold.GetProperties().volume > 0.0) {
output.push_back(outputManifold);
}
verts.clear();
}
} while (vl.inc());
return output;
}
std::vector<Manifold> Manifold::Fracture(
const std::vector<glm::vec3>& pts,
const std::vector<float>& weights) const {
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd like to avoid having two APIs for different precisions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it okay to expect everything to converge to double only with #659 ?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO yes.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes

std::vector<glm::dvec3> highpVerts(pts.size());
std::vector<double> 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> Manifold::ConvexDecomposition() const {
// Step 1. Get a list of all unique triangle faces with at least one reflex
// edge
std::unordered_set<int> 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::vec3 tangent =
glm::cross(impl.faceNormal_[faceA],
impl.vertPos_[impl.halfedge_[i].endVert] -
impl.vertPos_[impl.halfedge_[i].startVert]);
float tangentProjection = glm::dot(impl.faceNormal_[faceB], tangent);
// If we've found a pair of reflex triangles, add them to the set
if (tangentProjection > 1e-6) {
uniqueReflexFaceSet.insert(faceA);
uniqueReflexFaceSet.insert(faceB);
}
}
std::vector<int> uniqueFaces; // Copy to a vector for indexed access
uniqueFaces.insert(uniqueFaces.end(), uniqueReflexFaceSet.begin(),
uniqueReflexFaceSet.end());

// Step 2. Calculate the Circumcircles (centers + radii) of these triangles
std::vector<glm::dvec3> circumcenters(uniqueFaces.size());
std::vector<double> circumradii(uniqueFaces.size());
Vec<glm::dvec3> joggledVerts(impl.vertPos_.size());
for (size_t i = 0; i < impl.vertPos_.size(); i++) {
joggledVerts[i] = impl.vertPos_[i];
}
for (size_t i = 0; i < uniqueFaces.size(); i++) {
glm::dvec4 circumcircle = impl.Circumcircle(joggledVerts, uniqueFaces[i]);
circumcenters[i] =
glm::dvec3(circumcircle.x, circumcircle.y, circumcircle.z);
circumradii[i] = circumcircle.w;
}

// Step 3. If any two circumcenters are identical, joggle one of the triangle
// vertices, TODO: and store in a hashmap
std::mt19937 mt(1337);
std::uniform_real_distribution<double> dist(0.0, 1.0);
double randOffset = 0.0;
for (size_t i = 0; i < circumcenters.size() - 1; i++) {
for (size_t j = i + 1; j < circumcenters.size(); j++) {
if (glm::distance(circumcenters[i], circumcenters[j]) < 1e-6) {
joggledVerts[impl.halfedge_[(uniqueFaces[i] * 3) + 0].startVert] +=
glm::dvec3(dist(mt), dist(mt), dist(mt)) * 1e-11;
}
}
}

// Step 4. Recalculate the circumcenters
for (size_t i = 0; i < uniqueFaces.size(); i++) {
glm::dvec4 circumcircle = impl.Circumcircle(joggledVerts, uniqueFaces[i]);
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we joggle the circumcircles instead of the verts? Overriding vertPos_ in an Impl function is pretty strange.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My ultimate consideration is to record which verts were joggled before computing the voronoi diagram, and then to unjoggle them afterward so that they occupy their original exact coordinates again...

I'm considering alternative strategies, like applying an invertible spatial distortion to all the input vertices, and then backwards to all of the output voronoi diagram vertices.... just have to find the right strategy that preserves precision...

circumcenters[i] =
glm::dvec3(circumcircle.x, circumcircle.y, circumcircle.z);
circumradii[i] = circumcircle.w;
}

// Step 5. Calculate the Voronoi Fracturing
std::vector<Manifold> output = Fracture(circumcenters, circumradii);

// TODO:
// Step 6. Unjoggle the voronoi region vertices
// Step 7. Hull and Intersect with the original Manifold

return output;
}
} // namespace manifold
6 changes: 6 additions & 0 deletions src/third_party/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
1 change: 1 addition & 0 deletions src/third_party/voro
Submodule voro added at 56d619
18 changes: 18 additions & 0 deletions test/manifold_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -767,3 +767,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<Manifold> convexParts = nonConvex.ConvexDecomposition();

EXPECT_EQ(convexParts.size(), 224);

float originalVolume = nonConvex.GetProperties().volume;
float convex_volume = 0.0;
Manifold manifold_union = Manifold::Manifold();
for (Manifold cur_manifold : convexParts) {
manifold_union += cur_manifold.Hull();
}
float union_volume = manifold_union.GetProperties().volume;
EXPECT_NEAR(originalVolume, union_volume, 1e-6);
}
Loading