Manifold Performance

[elalish]

As a partner to #340, I'm starting this discussion as a place to share real-world performance insights on Manifold. I feel a discussion is a better form since performance varies widely based on the chosen benchmark, the parallel architecture, and what library it's compared against, not to mention the fact that our performance is improving over time.

To start, there's the first performance analysis I did, comparing Manifold to baseline OpenSCAD (CGAL), though Manifold has been optimized significantly since then (and likely OpenSCAD has as well). Overall takeaway was Manifold being 100 - 1,000 times faster.

Please post links to any performance analysis you've done and give a top-level summary here. I'm hoping we can build a more complete picture and also get a sense of what optimizations we should focus on going forward.

[elalish]

@ochafik made this excellent performance comparison of OpenSCAD running Manifold compared to CGAL's fast-csg, finding a 5-30x improvement. And fast-csg is itself a vast improvement over CGAL's Nef routines (baseline OpenSCAD), showing its own 30-150x improvement. This validates that the 1,000x improvement shown above was not a unique case.

[elalish]

There is also this lovely analysis comparing Manifold to JSCAD. Some quotes:

I then spent the rest of the day trying to figure out if it was actually true. I've been through it thoroughly and it appears that Manifold's booleans are 3 ORDERS OF MAGNITUDE (1000x) faster than the BSP-based booleans in JSCAD/CSharpCAD.

Using separate programs, I checked the output from Manifold and it appears correct. Indeed, in the case "union overlap", CSharpCAD was producing output that was not manifold. Manifold handled it correctly.

[pca006132]

Some of my previous Issues/PRs regarding performance:

• Performance tracking #105 and Optimizations #106
  Basically about eliminating unnecessary copying and special case faces with 3/4 vertices for faster processing.
• universal vector #121
  Switch to using unified memory instead of manually managing host and device vectors. This is faster because it can perform on-demand paging and doesn't require explicit copying between host and device vectors. This also allows new GPUs to over-commit memory and not fail when the mesh requires more memory than the VRAM available. This PR found a performance issue in thrust implementation of universal vector, and we resolved it by writing our own universal vector.
• dynamic backend implementation #126
  Enables backend selection based on the workload size, as multicore/GPU acceleration have certain overhead, and may not be beneficial if the workload is too small. This is related to Improve execution policy heuristics #380 and can be further improved.
• feature request: lazy boolean operations #114
  This is the current csg_tree module, which implements lazy boolean operation, tree rewriting and optimized batch boolean by changing evaluation order to minimize the mesh size. Basically stealing idea from Exploration of even lazier unions (bubbling through modules & transforms) for large speed ups openscad/openscad#3637 :P
• optimize collision detection #235
  Does more work by doing collision detection twice, and avoids atomic operations. This turns out to be faster than doing atomic operations.
• cache intermediate results (#199) #206 and Cache intermediate results automatically #199
  Fixes performance regression by adding a cache. Also related to Simplify the CsgOpNodes as we build them, rather than in GetChildren #368 which simplified the flattening logic and avoids stack overflow.

[pca006132]

Ideas about further optimization:

1. Parallelize CollapseEdge. This operation is the rare part that is not parallelized and takes a long time for large meshes. I think the operation should be local, so edges where the closed neighborhood of their end points are disjoint should be able to operate indepently? They should not touch even after collapsing the edges. (all edges may form a complicated graph which can be broken into connected components and be processed independently)
2. Use uint64_t for things like SparseIndices and morton code + id sorting. The uint64_t can work as a pair of uint32_t, and sorting performance with uint64_t can be 3 times faster than sorting zipped iterator.
3. Better parallelization for csg_tree, e.g. task based.

[pca006132]

The problem with picking random halfedges is that you still need a way to know which vertices (edges) are processed by the other threads now, and this requires some form of independent set.

I tried implementing a parallel connected component algorithm using union find, but the cost of calculating this is still too high (100ms when you have tens of millions of vertices + some book keeping to list out the connected components). Also, the problem is that for Sponge4, there are so many flagged edges that we only have 1 connected component even after filtering out those non-flagged edges. I tried splitting the edges into different groups (mod 5), and process 1 group each time, but this is also too slow.

[pca006132]

I also went ahead and implemented SparseIndices using 1 simple int64_t array (and modified a bunch of things to avoid using SwapPQ).

The sorting performance is improved at the expense of worse cache performance for some operations in boolean3.cpp, so the result is mixed when tbb is enabled. Some of the kernels only need p/q array, but because of how the array is now encoded, they need to access twice the memory.

For single threaded performance, we got around 10% performance improvement because we are not memory bound when using a single thread. Not sure if I should open a PR for this as the result is not as good as I hoped.

[elalish]

For this kind of thing, I would look at it from a code complexity standpoint. If the perf gain is not huge, I'd rather have code that's easier to read and maintain. What do you think?

[pca006132]

Regarding the parallelization of CollapseEdge, I am thinking about a spatial approach, similar to constructing a kd tree.
The idea is that we try to split vertices with their neighbors, only perform collapse on edges where the closed neighborhood is within the bounding box. When the two children are finished, we can merge the two bounding boxes (separated by a plane) and continue working on the remaining edges, and so on.

So the algorithm will be performed in two phases: The first phase is parallel and recursively splitting the vertices into a binary tree, where the pivot is the median of a fixed-size random set of vertices. This is kind of like quicksort? The second step is to perform the computation I mentioned above on each leave, and then merge two children leaves together and form a new leave, kind of like the merge phase of merge sort.

Implementation wise, this should be pretty efficient as we only need two buffers with size n, the binary tree is laid flat in memory, very similar to sorting. And we will also need to modify some parts of the code to avoid changing global states, e.g. push_back in FormLoop. Instead, we should store them into a backing store and increment an atomic counter for example.
I think this should be faster than sequentially collapsing all edges when we have a lot of flagged edges, e.g. > 100k or something. I haven't yet implemented this, but I tried running sorting algorithm and compare the time, the time of the sorting algorithm is < 1/20 of the overall time required for CollapseEdge if we have many edges to collapse.
We can figure out the number of flagged edges really quickly by changing to using tbb, as tbb allows using thread local constructs which is very useful for counting stuff.

[elalish]

Agreed, that sounds like a good approach. Our Collider in fact already does that - it creates a BVH to do spatial queries against. I bet you could add a method there to get the buckets you want. Let me know if the code is unclear - I wrote that a long time ago and it's probably a little more clever than it ought to be.

[pca006132]

Just to report on the hashtable changes: I tried implementing SparseMap with Hashtable for boolean3, and sorts are mostly removed (except for the sort for winding numbers). The problem is that Hashtable is really slow comparing with our previous array approach (basically double the time). I tried various approaches including

1. Resize the hashmap when the FIlledFraction is too low (< 1/8)
2. Optimize the size of hashmap. I found that changing the size from 2^n to 2^n-1 can reduce collision slightly, which results in large improvement in running time (>10%)
3. Relaxed atomic add memory order to relaxed. (not very helpful on x86_64 CPUs because atomic operations here must be sequential consistent)

While the time for Intersect12 is significantly faster, the time required for hashtable insert is just too slow. There are several reasons:

1. Atomic operations are just too slow. On x86_64, lock inc (atomic increment by 1) has a latency of (17-42), while the latency for typical instruction is <= 5. In addition to this, atomic operations will block (prevent reordering) a lot of things due to memory order requirements.
2. Memory locality does not exist for hash tables, leading to poor cache performance. (I tried using the identity function as hash function, the result is a lot worse due to high collision)
3. Significant false sharing issue. Basically, each cache line contains a set of values (usually 32 bytes), and CPU cores determine whether or not they are exclusive to itself or shared to other CPUs. When an atomic operation is done, the CPU has to do a lot more work when the data is shared across cores. For hashtable, we basically insert the entries into the hashtable in a pretty random way, so when we perform this in parallel we are having a huge false sharing issue.

And I think these issues are quite fundamental for hashtable, so I don't think there is a way out of this. We probably need to stay with the existing method.

[elalish]

Wow, thanks for the detailed analysis! Goes to show I still have some work to do in my performance intuition, though I do feel good that my early solution is holding up at least. Should we just leave it as is for now? I'm a bit concerned about thrust::unique, but perhaps it's better to wait a year until the pstl is available everywhere. What do you think?

[pca006132]

I feel that we can change unique to std sequential implementation. I guessed we are not gaining much from a parallel unique, and benchmark shows no measurable difference.

In fact, our performance mostly depends on parallel for and sort. Everything else is not too significant (will not double the running time)

[elalish]

That sounds ideal. Thanks!

[briansturgill]

I've been asking performance questions that you may find a bit odd... like why would I think simple ops take a while.
Here is something I commonly did in OpenSCAD... basically it's a way to fillet.
When I run the OpenSCAD below, it takes 18.75 seconds to render in CGAL.
When I run the equivalent code in Python using Manifold, it takes 0.228 seconds. Not a typo.

In OpenSCAD:

$fn = 100;

module rounded_rect(r, dim) {
    hull() {
        translate([r,r]) circle(r);
        translate([dim.x-r,r]) circle(r);
        translate([r,dim.y-r,r]) circle(r);
        translate([dim.x-r,dim.y-r]) circle(r);
    }
}

module rounded_slice(r, dim, h) {
    linear_extrude(h) rounded_rect(r, dim);
}

module rounded_slices(r, dim) {
    h = 0.1;
    union() {
        for(iota = [0.0: 0.1: 2.0]) {
            translate([iota, iota, iota]) rounded_slice(r-iota,[dim.x-2*iota, dim.y-2*iota], h);
        }
    }
}


rounded_slices(5, [10, 20]);

In Python:

import time
import manifold3d as m

def rounded_rect(r, dim):
    x, y = dim
    circ = m.CrossSection.circle(r, 100)
    return m.CrossSection.batch_hull(
        [
            circ.translate((0, 0)),
            circ.translate((x-r, 0)),
            circ.translate((0, y-r)),
            circ.translate((x-r, y-r)),
        ]
    )

def rounded_slice(r, dim, h):
    return m.Manifold.extrude(rounded_rect(r, dim), h)

def rounded_slices(r, dim):
    x, y = dim
    h = 0.1
    l=[]
    iota = 0.0
    while iota < 2.09:
        l.append(
            rounded_slice(r-iota, (x-2*iota, y-2*iota), h).
            translate([iota, iota, iota])
            )
        iota += 0.1

    return m.Manifold.batch_boolean(l, m.OpType.Add)

ts = time.perf_counter()
rs =rounded_slices(5, (10, 20))
rs.num_vert() # Make sure all is evaluated.
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "Created round_slices created.")

import piecad
piecad.view(piecad.Obj3d(rs))

[briansturgill]

@pca006132 I'm really glad to see that, but I really want to be doing an instructional/maker oriented CAD package in Python.
When I was making CSharpCAD, I had many ideas. I want to bring those out. Literally all CAD packages have too
much jargon and and oddities like making you worry about winding order. I have ideas for making that easier.
Alas with CSharpCAD, I was unable to solve the 3D booleans issue. I had a series of misadventures involving health,
the death of my father and moving that made me abandon CSharpCAD for a time. When I finally got back, I found
Manifold had added Clipper2. I also found I didn't like actually making models in C#... I always ended up using
the Python front-end. Thus I'm here.

The two most widely known programming languages are JavaScript and Python.
Python is the most widely used in instruction though exposure to JS is usually also given.
Amongst makers, Python is hands down better known.
Makers generally just want to get something done. It may have been months since the last time they did something
with CAD. Again, Python plus a simple API is the answer.

OpenSCAD is at least somewhat simple. But its functional nature makes many things far harder than they need to be.
For example, for many years you had to do unbelievable tricks to output multiple single objects from the same script.
Even now, it's really strange what you have to do. It is very hard to debug OpenSCAD. I have something I call
cadviewer that is a separate program that gets (via an HTTP POST) a 3d object you want to see. You send an object
with the view function. Each view is captured and retained by cadviewer. You can uses arrow keys to switch between
objects. If you need to capture each return of a functions, you can put something like return view(obj), view
returns it's first argument. Combine that with the great Python debugging experience many IDEs have and
OpenSCAD loses hands down. IDEs are available for all levels of sophistication. OpenSCAD tying you into a
single development environment is bad. Most CAD users do CAD rarely, they need to be able to quickly come up
to speed each time and use of their favorite environment is a big help here.

Note, cadviewer can easily be used from JSCAD or Manifold C++... anywhere that can do an HTTP POST of a 3D object.
I handle 2D objects by having view extrude it thinly.

[elalish]

Yeah, I agree, hence why we developed ManifoldCAD.org. I even made a Python colab example to show how to use Python to script objects in a shareable web interface. The sky's the limit, whether you want to collaborate/extend some of these projects or make your own front-end: https://github.com/elalish/manifold?tab=readme-ov-file#manifold-frontend-sandboxes

[briansturgill]

I tried running the Colab notebook: Both at google and on my local system, I get:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: object of type 'int' has no len()

The above exception was the direct cause of the following exception:

SystemError                               Traceback (most recent call last)
Cell In[2], line 89
     87 # Main
     88 t0 = time()
---> 89 manifold = halfSponge(3)
     90 t1 = time()
     91 print(f"Took {(t1-t0)*1000:.1f}ms")

Cell In[2], line 84, in halfSponge(n)
     81 result -= hole.rotate([90, 0, 0])
     82 result -= hole.rotate([0, 90, 0])
---> 84 return result.trim_by_plane([1, 1, 1], 0).set_properties(4, posColors).scale(100)

SystemError: <nanobind.nb_method object at 0x7fe8278aba00> returned a result with an exception set

[briansturgill]

So I removed halfSponge from the picture (just made a cube).
That's really cool I didn't realize that trimesh had that level of Jupyter support.

[elalish]

Oh man, you're right! Looks like the python sponge example got updated for some API fix and I didn't update the colab. Should be fixed now.

[briansturgill]

I ran several timing tests to determine several things.

1. Is it practical to create heavily used objects from Python?
   Answer: Yes, the overhead from Python needing to pass points is minimal]
2. Explore issues relating to a problem I had making rounded project boxes in OpenSCAD .(Manifold Performance #383 (comment))
   Answer: Manifold wins here by more than 2 orders of magnitude. Even 2D hulls are practial.
3. 3D Hulls... like any place I've seen, not stellar in performance, but definitely acceptable.
4. Speed of Unions and BatchUnions
   Answer: 2d fine, especially as you rarely do very many booleans there. 3d stellar.
5. Mistake in earlier listing... I finally got the 2d and 3d objects to finalize correctly.
   Some corrections of earlier statements... as numbers have significantly changed.
   a) While 2D batch_boolean is still significantly faster than individual unions, 3d batch_boolean isn't.
   b) Manifold::Cube is slower than using a Python function that uses square and extrude. Manifold::Cube is slow. #710
   [Sorry, but I do find that funny! :-)]

2D timing tests

 0.481 1,000,000 creations of CS with square member function
 1.953 1,000,000 creations of CS with points from square
 1.136 100,000 creations of roundrect using hull and 4 circles
 0.002 Creation of 1,000 randomly placed roundrect
 1.130 1000 unions of 1000 randomly placed roundrect
 0.522 batch_boolean union of 1,000 randomly placed roundrect

import time
import manifold3d as m
from random import random

ts = time.perf_counter()
for i in range(1000000):
    c = m.CrossSection.square((10.0, 10.0))
    c.num_vert()
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "1,000,000 creations of CS with square member function")

p = c.to_polygons()
ts = time.perf_counter()
for i in range(1000000):
    c = m.CrossSection(p)
    c.num_vert()
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "1,000,000 creations of CS with points from square")


r = 5
s = 36
ts = time.perf_counter()
for i in range(100000):
    circ = m.CrossSection.circle(r, s)
    c = m.CrossSection.batch_hull(
        [
            circ,
            circ.translate((10, 0)),
            circ.translate((0, 10)),
            circ.translate((10, 10)),
        ]
    )
    c.num_vert()
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "100,000 creations of roundrect using hull and 4 circles")

l = []
ts = time.perf_counter()
for i in range(1000):
    c1 = c.translate((random() * 200, random() * 200))
    c1.num_vert()
    l.append(c1)
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "Creation of 1,000 randomly placed roundrect")

c1 = m.CrossSection()
ts = time.perf_counter()
for i in range(1000):
    c1 = m.CrossSection.__add__(c1, l[i])
    c1.num_vert()
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "1000 unions of 1000 randomly placed roundrect")

ts = time.perf_counter()
c = m.CrossSection.batch_boolean(l, m.OpType.Add)
c.num_vert()
te = time.perf_counter()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "batch_boolean union of 1,000 randomly placed roundrect")

3D timing tests

 1.291 100,000 creations of M with cube member function
 1.298 100,000 creations of M with mesh from cube
 0.719 100,000 creations of M with python function
 3.393 1,000 creations of roundedbox using hull and 8 spheres
 0.072 Creation of 1000 randomly placed roundedbox
 10.466 1000 unions of 1000 randomly placed roundedbox
 9.868 batch_boolean union of 1000 randomly placed roundedbox

import time
import manifold3d as m
from random import random
import numpy as np

ts = time.time()
for i in range(100000):
    c = m.Manifold.cube((10.0, 10.0, 10.0))
    c.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "100,000 creations of M with cube member function")

mesh = c.to_mesh()
tv=np.array(mesh.tri_verts.tolist().copy(), np.int32)
vp=np.array(mesh.vert_properties.tolist().copy(), np.float32)
ts = time.time()
for i in range(100000):
    mesh = m.Mesh(tri_verts=tv, vert_properties=vp)
    c = m.Manifold(mesh)
    c.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "100,000 creations of M with mesh from cube")


def cube(dim):
    x, y, z = dim
    c = m.CrossSection.square([x, y])
    return m.Manifold.extrude(c, z)


ts = time.time()
for i in range(100000):
    c = cube((10.0, 10.0, 10.0))
    c.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "100,000 creations of M with python function")


r = 5
s = 36
ts = time.time()
for i in range(1000):
    sph = m.Manifold.sphere(r, s)
    c = m.Manifold.batch_hull(
        [
            sph,
            sph.translate((10, 0, 0)),
            sph.translate((0, 10, 0)),
            sph.translate((10, 10, 0)),
            sph.translate((0, 0, 10)),
            sph.translate((10, 0, 10)),
            sph.translate((0, 10, 10)),
            sph.translate((10, 10, 10)),
        ]
    )
    c.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", "1,000 creations of roundedbox using hull and 8 spheres")

n = 1000
l = []
ts = time.time()
for i in range(n):
    c1 = c.translate((random() * 200, random() * 200, random() * 200))
    c1.num_vert()
    l.append(c1)
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", f"Creation of {n} randomly placed roundedbox")

c1 = m.Manifold()
ts = time.time()
for i in range(n):
    c1 = m.Manifold.__add__(c1, l[i])
c1.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", f"{n} unions of {n} randomly placed roundedbox")

ts = time.time()
c = m.Manifold.batch_boolean(l, m.OpType.Add)
c.num_vert()
te = time.time()
tdiff = te - ts
print(f"{tdiff: 5.3f}", f"batch_boolean union of {n} randomly placed roundedbox")

[pca006132]

For performance issues, you probably want to look at https://github.com/elalish/manifold/wiki/Performance-Considerations. Feel free to let me know if you want to add something there. Basically, python is not openscad, you have to do caching manually. The more you cache, the better the performance (in general).

For 3D hull, we are using an external library. It seems to be the best library we can find, but there are also some issues such as it may print some error messages to stderr, and it is not parallelized. We may be able to improve its performance a bit, but currently no one is working on that. Maybe it can be a good GSoC project as well.

[hrgdavor]

I have no real world performance numbers here, but just and idea.

Would it be beneficial and generally useful for users using extrude to provide the 2d shape, and already triangulated mesh for it instead of relying on manifold to triangulate it ?

[pca006132]

Not really, we have a lot more triangulation going on in boolean operations, extrude is rarely the slowest part.

[OLDMAN75914091]

Is there a generally agreed-upon definition of Performance? To most it means execution speed. That's fine, but speed at the expense of robustness/reliability is no gain. Maybe I missed it in my readings, but there have been precious few discussions on reliability/correctness of Manifold's Boolean output. SOTA Boolean libraries with a legacy of many decades are still struggling with this general issue. ParaSolids from Siemens is one. PTC and Dassault Systemes have their own. But, I believe the issue remains. Which explains my curiosity on how/if Manifold addresses this.

[elalish]

There's a polygon corpus in our tests you can use. Lots of ugly things from real Boolean intersections in there.

Take a close look at a CAD fillet at the intersection of two curved NURBS surfaces and let me know how precise it is. :)

I think our curved triangles have some potential as lossy compression by fitting to a bunch of verts. We'll see.

[OLDMAN75914091]

where do i find this polygon corpus ?

[pca006132]

https://github.com/elalish/manifold/tree/master/test/polygons

[OLDMAN75914091]

" I think our curved triangles have some potential as lossy compression by fitting to a bunch of verts. We'll see. "

Yes, that could well be.
But since it is going to be lossy, it might not be suitable for CAD--->CAE--->CAM

[elalish]

True, for CAD/CAM I think it might be more useful as a way to create manifold fillets from chamfers, though I still need to build a chamfer function.

[OLDMAN75914091]

You mentioned "symbolic perturbation" (of vertices, I suppose).
I understand the idea behind perturbation of vertices (usually for corner cases)
But what does symbolic mean?

[elalish]

It means they don't actually move, hence why we avoid corner cases entirely. In our case, instead of e.g. a vert is on one side of a plane or the other, or exactly on (<, >, =), we use only (<=, >) - the equal case is symbolically perturbed to one side, so nothing is ever "on".

[OLDMAN75914091]

Stumbled into this libary recently: https://meshlib.meshinspector.com/

Have you folks checked out their 3D Booleans?

[elalish]

No, I hadn't heard of them, thanks for sharing! I'm pleased they listed us as their top competitor for speed. 😁 I would love to see how they do on our test cases; based on what I see in their code, I would guess they do not have the same manifoldness guarantee as us.