Add MinkowskiSum and MinkowskiDifference #666
Conversation
|
The non-convex/non-convex test appears to be failing on some of the GH Actions test runners due to triangulations not coming out CCW... It passes on Linux (and my personal Windows machine), but fails on Github's Mac and Windows Runners. I wonder if these will pass on the We can also say that Non-Convex/Non-Convex is not supported, but that would make me sad... |
|
It seems that there is something weird going on:
The epsilon is small (7e-6), no way this is a valid polygon, and I don't think it is caused by precision issue. |
How did you generate this image? Since this implementation is just unioning hulls of triangle pairs, I believe it's acting like a fuzzer for the union operation, where every hull is the worst case of "just barely touching" the other hulls. Whatever assumptions about floating point Manifold makes seem to hold true on some platforms, but not so much on others... One mysterious thing to note is that the number of triangles output by both the sequential and parallel unions is not constant between runs; I think there might be an indeterminism sneaking into the library somewhere... I'll force the CI to run again to see if different machines pass and fail... |
|
We have non-determinism by design. We can force determinism by setting The plot: import numpy as np
import matplotlib.pyplot as plt
array = np.array
show = lambda x: plt.plot(*np.hsplit(np.append(x, x[0, :].reshape((1, 2)), axis=0), 2), marker='.')
show(array([
[0.0500000268, 0.435410291],
[-0.00156550109, 0.413770884],
[0.0366442874, 0.405546039],
[0.0366442874, 0.405546039],
[0.0500000268, 0.402671158],
]))
show(array([
[-0.0454634838, 0.395349145],
[-0.0454634838, 0.395349145],
[0.0500000268, 0.435410291],
[-0.0131011149, 0.397831321],
[0.0500000268, 0.402671158],
[0.0366442874, 0.405546039],
[0.0366442874, 0.405546039],
[-0.00156550109, 0.413770884],
]))
plt.show() |
|
It is interesting to note that different jobs are failing now, which suggests it might be a quirk of the order of operations... I don't suppose Manifold can catch itself as soon as it generates an inconsistent triangulation, and dump the two manifolds that it tried to union just before doing so? There's also a chance that the test is just too expensive for the runner, and it's inconsistently running out of memory 🤔 |
|
OOM: I don't think that will cause this behavior.
Actually we can do that with some hacking, need to modify |
|
Thanks for this! I've been toying with a similar idea but just for offsetting (minkowski with a sphere). I'll put up a PR soon so we can compare and contrast. Anyway, I ran into I believe the same trouble as you and came to the same conclusion - that this is a good fuzz test of our Booleans. I have a hunch that I see two Boolean bugs I may be able to fix - one regarding our symbolic perturbation and another that's creating those self-intersecting polygons. |
|
By the way, there’s a neat trick with Minkowski Sums: “Linear” Shape Interpolation
It’s the old lerp formula This method even works on shapes with different Genuses. Not sure what it can be used for in modeling just yet 😅 |
|
Okay, that is pretty cool - Minkowski is growing on me. |
|
GLM is header-only - I'd be a little surprised if it produced its own lib? By main system unit tests, do you just mean the Minkowski tests you added? That means you're producing self-intersecting polygons (or other degenerate geometry the triangulator is having trouble with). Best would be if you can reduce them to a minimal Boolean op that fails with the fewest vertices; then we can debug. |
Perhaps there's an issue with the default CMake config? I notice that I can't get it to set Perhaps I could see a working Windows CMake config as an example? Maybe we can fix the first-time automatic setup...
Yeah; just the new tests. I haven't been able to look at the output yet, but the op is mostly just the union of a bunch of hulls, so it shouldn't be too heinous... Once I've got it running locally, it'll be easier to isolate. |
|
Here's what the CI is doing, which seems to work: https://github.com/elalish/manifold/blob/master/.github/workflows/manifold.yml#L180 |
Son of a gun; I notice it's on in the default config; I wonder if it's doing good work there, or if we should make it a Windows-specific toggle... |
|
Well... it doesn't seem like the switch to But the simplification does resolve a lot of the crunchiness, leading to much cleaner topology in the end 😄 cube = manifold3d.Manifold.cube([1.0, 1.0, 1.0], True)
fun_shape = cube - cube.translate([0.5, 0.5, 0.5])
small_sphere : manifold3d.Manifold = manifold3d.Manifold.sphere(0.1, 25)
small_sphere2 : manifold3d.Manifold = manifold3d.Manifold.sphere(0.2, 25)
general_offset = fun_shape .minkowski_sum (small_sphere)
general_offset = general_offset.minkowski_difference(small_sphere2)
general_offset = general_offset.minkowski_sum (small_sphere)Old: New: Old Rotated: New Rotated: I did disable parallelization here because it seems like Thrust was removed and I'm not sure how to replace it... |
|
I'm not sure what you mean by genus non-determinism in the CI - I'm seeing mostly triangulations failures and seg faults. Indeed, Thrust is out, but we already had wrapper functions for things like |
|
In general, it is really hard to make things deterministic while having parallelization with floating point operations. I think we had some basic supports, but it requires enabling some option, and it will make things slower. |
I see; thank you for pointing that out to me. I parallelized the Convex+NonConvex case according to this recommendation. Are there any instances parallelizing a concurrent append buffer? The NonConvex+NonConvex case could probably benefit a lot from parallelization, but I'm not sure how to do it with the existing primitives... std::vector + mutex maybe?
I think I had the CI test misconfigured; the genus + surface area appears to be stable now on Linux and my home Windows machine; just the CCW issues on Mac, CI Windows, and TBB Linux 🤔 Quite suspicious it's only those CI systems...
I've enabled the option for the CI test, for what it's worth. Hopefully I don't actually have determinism issues in this regard. 😅 |
|
I am quite interested in the failing cases with non-CCW triangles. I think those cases reveal some issues in our boolean algorithm? |
|
E.g. this can't be a valid polygon with precision 1e-12:
|
|
Maybe we should add some code to dump the offending meshes out when the triangulation fails? |
This reverts commit d1f35e2.
|
Ah, interesting: I can reproduce the CI error if I enable
I copied the code for it from other tests and tried, but I haven't figured out how to get assimp and CMake to play nicely with each other yet 😅 I can also replicate the tests in Python... |
|
Here is a Python reproduction of the test, with the mesh outputs that are supposedly causing these triangulation errors: My bet is that it's partially overlapping, non-axis-aligned, coplanar faces, like here in Though, it's interesting that Prusa Slicer doesn't detect anything wrong with these meshes... It's also interesting to me that as_original doesn't seem to be applying the simplification to this one... I think because it's called before it is lazy constructed? Would it make sense to force minkowski mesh construction early so that simplification can be called on it? |
|
OK, I have reduced the failure case to TEST(BooleanComplex, TriangulationFail) {
// clang-format off
MeshGL64 m;
m.precision = 2e-12;
m.vertProperties = {
0.75, -1, -1, //
0.75, -1, -0.25, //
0.75, -0.666666666666667, 0.0833333333333333, //
0.75, -0.25, -1, //
0.75, 0.0833333333333333, -0.666666666666667, //
0.833333333333333, 0.5, -0.166666666666667, //
0.75, 0.0416666716337204, 0.0416666716337204, //
0.75, 0.5, 0.25, //
0.75, 0.25, 0.5, //
1.25, -1, -1, //
1.25, -1, -0.75, //
1.25, -1, -0.75, //
1.25, -0.75, -1, //
1, -0.5, 0, //
1, -0.25, -0.25, //
1, 0, -0.5, //
1.25, 0.5, 0.25, //
1.25, 0.5, 0.75};
m.triVerts = {
1, 6, 0, //
3, 0, 6, //
0, 10, 1, //
0, 3, 12, //
1, 13, 2, //
6, 1, 2, //
13, 8, 2, //
15, 3, 4, //
4, 3, 6, //
15, 4, 7, //
5, 15, 7, //
6, 2, 8, //
6, 14, 7, //
8, 14, 6, //
7, 4, 6, //
13, 14, 8, //
5, 7, 16, //
0, 9, 10, //
9, 0, 12, //
1, 10, 13, //
12, 3, 15, //
13, 10, 14, //
11, 7, 12, //
12, 14, 10, //
10, 9, 12, //
17, 11, 12, //
14, 12, 7, //
15, 5, 16, //
16, 12, 15, //
7, 11, 17, //
7, 17, 16, //
12, 16, 17};
MeshGL64 m2;
m2.precision = 2e-12;
m2.vertProperties = {
0.75, 0, 0.75, //
0.75, 0.25, 0.5, //
1.25, -1, -0.75, //
1.25, -0.75, -1, //
1, 0.5, 0.5, //
1.25, 0.5, 0.25, //
1.25, 0.5, 0.75};
m2.triVerts = {
4, 1, 0, //
3, 0, 1, //
0, 3, 2, //
3, 6, 2, //
0, 2, 6, //
3, 1, 5, //
1, 4, 5, //
4, 0, 6, //
3, 5, 6, //
4, 6, 5};
// clang-format on
ExportMesh("trifail.glb", (Manifold(m) + Manifold(m2)).GetMeshGL(), {});
} |
|
I'm having a little trouble following everything that's going on here. Yes, without |
|
Unfortunately, I don't get a repro with that test - it seems to be working on mac. Might be good to put it in its own PR so we can see what happens on the CI. I'll just do that now. |
|
Okay, got the repro to work, but that's already too late - the inputs themselves are already self-intersecting. We'll need to find a way to determine earlier when something goes wrong. Might be a good time to write an approximate self-intersection check function. |
|
I see, I guess the simplest way would be to convert the mesh into CGAL mesh, and it will report an error when it finds self-intersections? The problem is that it will treat epsilon-valid meshes as self-intersecting, not sure if that is a problem here. |
|
Probably worth a try at least - we're just debugging. Ideally we shouldn't even be creating marginal geometry. |
|
Patch for debugging: diff --git a/src/manifold/CMakeLists.txt b/src/manifold/CMakeLists.txt
index 00136b5d..669188a1 100644
--- a/src/manifold/CMakeLists.txt
+++ b/src/manifold/CMakeLists.txt
@@ -20,9 +20,10 @@ add_library(${PROJECT_NAME} ${SOURCE_FILES})
target_include_directories(${PROJECT_NAME} PUBLIC
$<INSTALL_INTERFACE:include>
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>)
+find_package(CGAL REQUIRED COMPONENTS Core)
target_link_libraries(${PROJECT_NAME}
PUBLIC utilities
- PRIVATE $<BUILD_INTERFACE:collider> polygon ${MANIFOLD_INCLUDE}
+ PRIVATE $<BUILD_INTERFACE:collider> polygon CGAL::CGAL CGAL::CGAL_Core
)
target_compile_options(${PROJECT_NAME} PRIVATE ${MANIFOLD_FLAGS})
diff --git a/src/manifold/include/manifold/manifold.h b/src/manifold/include/manifold/manifold.h
index c59a9149..4c3e40a7 100644
--- a/src/manifold/include/manifold/manifold.h
+++ b/src/manifold/include/manifold/manifold.h
@@ -321,10 +321,10 @@ class Manifold {
///@}
struct Impl;
+ Manifold(std::shared_ptr<Impl> pImpl_);
private:
Manifold(std::shared_ptr<CsgNode> pNode_);
- Manifold(std::shared_ptr<Impl> pImpl_);
static Manifold Invalid();
mutable std::shared_ptr<CsgNode> pNode_;
diff --git a/src/manifold/src/csg_tree.cpp b/src/manifold/src/csg_tree.cpp
index f62ba5cf..1d9c6477 100644
--- a/src/manifold/src/csg_tree.cpp
+++ b/src/manifold/src/csg_tree.cpp
@@ -26,9 +26,58 @@
#include "impl.h"
#include "manifold/parallel.h"
#include "mesh_fixes.h"
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
constexpr int kParallelThreshold = 4096;
+using NT3 = CGAL::Gmpq;
+using CGAL_Kernel3 = CGAL::Cartesian<NT3>;
+using Polyhedron = CGAL::Polyhedron_3<CGAL_Kernel3>;
+
+class CGALPolyhedronBuilderFromManifold : public CGAL::Modifier_base<typename Polyhedron::HalfedgeDS>
+{
+ using HDS = typename Polyhedron::HalfedgeDS;
+ using CGAL_Polybuilder = CGAL::Polyhedron_incremental_builder_3<typename Polyhedron::HalfedgeDS>;
+public:
+ using CGALPoint = typename CGAL_Polybuilder::Point_3;
+
+ const manifold::MeshGL64& meshgl;
+ CGALPolyhedronBuilderFromManifold(const manifold::MeshGL64& mesh) : meshgl(mesh) { }
+
+ void operator()(HDS& hds) override {
+ CGAL_Polybuilder B(hds, true);
+
+ B.begin_surface(meshgl.NumVert(), meshgl.NumTri());
+ for (size_t vertid = 0; vertid < meshgl.NumVert(); vertid++) {
+ auto pt = meshgl.GetVertPos(vertid);
+ B.add_vertex(CGALPoint(pt.x, pt.y, pt.z));
+ }
+
+ for (size_t faceid = 0; faceid < meshgl.NumTri(); faceid++) {
+ const auto tv = meshgl.GetTriVerts(faceid);
+ B.begin_facet();
+ for (const int j : {0, 1, 2}) {
+ B.add_vertex_to_facet(tv[j]);
+ }
+ B.end_facet();
+ }
+ B.end_surface();
+ }
+};
+
+namespace {
+void manifoldToCGALSurfaceMesh(Manifold &manifold) {
+ auto p = std::make_shared<Polyhedron>();
+ auto meshgl = manifold.GetMeshGL64();
+ CGALPolyhedronBuilderFromManifold builder(meshgl);
+ p->delegate(builder);
+ printf("%d\n", p->is_closed());
+ CGAL::Nef_polyhedron_3<CGAL_Kernel3> x(*p);
+}
+}
+
namespace {
using namespace manifold;
struct Transform4x3 {
@@ -458,8 +507,17 @@ std::shared_ptr<Manifold::Impl> CsgOpNode::BatchBoolean(
if (results.size() == 1)
return std::make_shared<Manifold::Impl>(*results.front());
if (results.size() == 2) {
- Boolean3 boolean(*results[0], *results[1], operation);
- return std::make_shared<Manifold::Impl>(boolean.Result(operation));
+ try {
+ Boolean3 boolean(*results[0], *results[1], operation);
+ auto r = std::make_shared<Manifold::Impl>(boolean.Result(operation));
+ auto m = Manifold(r);
+ manifoldToCGALSurfaceMesh(m);
+ return r;
+ } catch (manifold::geometryErr err) {
+ results[0]->Dump();
+ results[1]->Dump();
+ throw err;
+ }
}
#if MANIFOLD_PAR == 'T' && __has_include(<tbb/tbb.h>)
if (!ManifoldParams().deterministic) {
@@ -479,11 +537,19 @@ std::shared_ptr<Manifold::Impl> CsgOpNode::BatchBoolean(
continue;
}
group.run([&, a, b]() {
- const Manifold::Impl *aImpl;
- const Manifold::Impl *bImpl;
- Boolean3 boolean(*getImplPtr(a), *getImplPtr(b), operation);
- queue.emplace(
- std::make_shared<Manifold::Impl>(boolean.Result(operation)));
+ const Manifold::Impl *impla = getImplPtr(a);
+ const Manifold::Impl *implb = getImplPtr(b);
+ try {
+ Boolean3 boolean(*impla, *implb, operation);
+ auto r = std::make_shared<Manifold::Impl>(boolean.Result(operation));
+ auto m = Manifold(r);
+ manifoldToCGALSurfaceMesh(m);
+ queue.emplace(r);
+ } catch (manifold::geometryErr err) {
+ impla->Dump();
+ implb->Dump();
+ throw err;
+ }
return group.run(process);
});
}
@@ -507,12 +573,19 @@ std::shared_ptr<Manifold::Impl> CsgOpNode::BatchBoolean(
auto b = std::move(results.back());
results.pop_back();
// boolean operation
- Boolean3 boolean(*a, *b, operation);
- if (results.size() == 0) {
- return std::make_shared<Manifold::Impl>(boolean.Result(operation));
+ try {
+ Boolean3 boolean(*a, *b, operation);
+ auto r = std::make_shared<Manifold::Impl>(boolean.Result(operation));
+ auto m = Manifold(r);
+ manifoldToCGALSurfaceMesh(m);
+ if (results.size() == 0)
+ return r;
+ results.push_back(r);
+ } catch (manifold::geometryErr err) {
+ a->Dump();
+ b->Dump();
+ throw err;
}
- results.push_back(
- std::make_shared<const Manifold::Impl>(boolean.Result(operation)));
std::push_heap(results.begin(), results.end(), cmpFn);
}
return std::make_shared<Manifold::Impl>(*results.front());
@@ -573,8 +646,10 @@ void CsgOpNode::BatchUnion() const {
tmp.push_back(
std::dynamic_pointer_cast<CsgLeafNode>(children_[start + j]));
}
- impls.push_back(
- std::make_shared<const Manifold::Impl>(CsgLeafNode::Compose(tmp)));
+ auto result = std::make_shared<Manifold::Impl>(CsgLeafNode::Compose(tmp));
+ auto m = Manifold(result);
+ manifoldToCGALSurfaceMesh(m);
+ impls.push_back(result);
}
}
diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp
index 590e5186..49555999 100644
--- a/src/manifold/src/impl.cpp
+++ b/src/manifold/src/impl.cpp
@@ -731,4 +731,32 @@ SparseIndices Manifold::Impl::VertexCollisionsZ(VecView<const vec3> vertsIn,
return collider_.Collisions<false, false>(vertsIn);
}
+void Manifold::Impl::Dump() const {
+#ifdef MANIFOLD_DEBUG
+ std::cout << std::setprecision(15);
+ std::cout << "MeshGL64 m;" << std::endl;
+ std::cout << "m.precision = " << precision_ << ";" << std::endl;
+ std::cout << "m.vertProperties = {" << std::endl;
+ bool first = true;
+ for (const auto& pos : vertPos_) {
+ if (!first) {
+ std::cout << ", //" << std::endl;
+ }
+ first = false;
+ std::cout << " " << pos.x << ", " << pos.y << ", " << pos.z;
+ }
+ std::cout << "};" << std::endl;
+ std::cout << "m.triVerts = {" << std::endl;
+ for (size_t i = 0; i < NumTri(); i++) {
+ if (i != 0) {
+ std::cout << ", //" << std::endl;
+ }
+ std::cout << " " << halfedge_[3 * i].startVert << ", "
+ << halfedge_[3 * i + 1].startVert << ", "
+ << halfedge_[3 * i + 2].startVert;
+ }
+ std::cout << "};" << std::endl;
+ std::cout << "----" << std::endl;
+#endif
+}
} // namespace manifold
diff --git a/src/manifold/src/impl.h b/src/manifold/src/impl.h
index 6404aa81..05c6a9ab 100644
--- a/src/manifold/src/impl.h
+++ b/src/manifold/src/impl.h
@@ -358,5 +358,8 @@ struct Manifold::Impl {
// quickhull.cpp
void Hull(VecView<vec3> vertPos);
+
+ // debug
+ void Dump() const;
};
} // namespace manifoldCGAL seems to be happy so far, and then we get invalid polygon for triangulation:
The mesh is pretty large, so I am not attaching the input here. I think we can trim down the mesh size by intersecting with the bounding box of the smaller one. |
|
The mesh size is the smallest when test with |
This is the minimal code to get the naive minkowksi sum working
with basic affordances for multithreading.This is essentially the other PR, but without the in-progress experimentation or the controversial
voro++dependency.The benefits of the Naive Technique are that:
The drawbacks to the Naive Technique are that:
There are two basic versions of the function, one with significant threading and the other without.I expect we'll consolidate to one or the other as their respective pros/cons make themselves apparent; the speed appears to be pretty comparable on my windows machine (but perhaps that's not true on Linux?).
I've decided the speed difference is too negligible (1-2% at best), and it's now a source of failing unit tests. So no more threading!
Tests and Bindings have been added 😎