Mesh Offsetting

[pca006132]

Mesh offsetting seems very useful for adding tolerance to parts, especially when combined with face selection.

~~I wonder if I can just offset each vertices according to their tangent? Although we may need to recompute the tangent if we only selected a subset of faces connecting to that vertex.~~ OK this will not work.

[pentacular]

You run the risk of self-intersection.

Also, offsetting by the tangent of the vertex won't offset the adjacent face by a fixed amount, which may give uneven tolerance across the face.

[pca006132]

So maybe we can compute the offset of each selected faces and add them to the vertices? Not sure about self intersection though.

[elalish]

Yeah, I was thinking about adding the vertex normal as a const input to the Warp function to allow this kind of thing, but indeed, it is very easy for a warping function to introduce self-intersections. The "right" way to offset a mesh is to do a Minkowski operation with a sphere (OpenSCAD supports this) but it's very slow and dramatically explodes the number of faces.

[thehans]

To avoid self intersections, you need to guarantee that the input is convex. So what is needed is a convex decomposition operation (then offset each convex piece, then union to get the final result).

I wonder if I can just offset each vertices according to their tangent? Although we may need to recompute the tangent if we only selected a subset of faces connecting to that vertex. OK this will not work.

I don't think a vertex has any specific "tangent", maybe you mean the vertex normal? Vertex normal makes a bit more sense as that is a concept from the field of graphics. But i'd consider it more like a hack/fudged values that looks ok to the eye, than something of rigorously defined geometry. Imagine translating each face by the offset, along its normal; then conceptually extend each face until it once again intersects at new edges/vertices. That's where you want the vertex to end up, but I'm not entirely sure if the straight average of face normals, applied to the vertex results in this transformation or if there is more to it.

[elalish]

The reason a Minkowski sum is the "right" way is that for a general vertex (not symmetric, more than three faces), when you offset the faces, they will no longer meet up at that same single vertex (in any position). The simplest alternative is just to offset the triangles exactly along their face normals, then stitch the edge and vertex gaps with new faces, but again it only works on convex polyhedra and it explodes the number of faces rapidly.

An alternative approach to adjusting tolerance is to design things parametrically in the first place, so they can simply be rapidly regenerated to adjusted dimensions. This is the approach I always took with OpenSCAD, and it'll work equally well with Manifold.

I tried for some time (years ago) to create a 3D straight skeleton algorithm, which would solve mesh offsetting, but it turned out to be more difficult than I imagined. I'm not tempted to spend much time on it anymore, but I agree it would be a nice algorithm to have.

[zalo]

The dumbest possible Minkowski algorithm I know (for a convex “tool” object) is to go through each face of one mesh, hull the other mesh to each of the vertices of each face, and then union them all together with the original mesh.

The hulls are all parallelizable, there are a ton of duplicate vertices, and it’s one giant union… I wonder how fast it would be?

EDIT: Seems like there might be some simple exact convex decomposition algorithms out there too, though the hulling step here might be slow… https://www.reddit.com/r/GraphicsProgramming/s/0jGSbAPnxF

[pca006132]

The hulls are all parallelizable, there are a ton of duplicate vertices, and it’s one giant union… I wonder how fast it would be?

I think it will not scale well for meshes with over a thousand faces...

Seems like there might be some simple exact convex decomposition algorithms out there too, though the hulling step here might be slow…

We want to avoid exact convex decomposition not because the algorithm itself is hard, but because exact convex decomposition will create many parts for very simple geometries, e.g. cube - sphere. In that case it is not better than the "dumbest possible" solution you said.

I actually implemented something that can deal with concave geometries, see https://github.com/elalish/manifold/issues/415#issuecomment-1725063678

The reason I was holding it back is because

1. The algorithm itself will create self-intersections when the curvature is too large / offset is too large that requires changing the topology globally. So I need another algorithm to compute some kind of decomposition that can avoid the self-intersection, and I found a flaw in the other algorithm I wrote (not hard to fix, just need to rewrite some code).
2. It is doing some approximation for concave surfaces. It makes me feel a bit uncomfortable...
3. We should round over those convex edges and vertices, similar to minkowski with a sphere with high $fn. But this is tricky and the paper I read just hand waived over this. I experimentally figured out a solution that seems to work. While the idea is simple, it is not very easy to implement, especially when I have some "clever" (complicated) index manipulation within the code...
4. It is unclear to me how I should clean this up and submit in a PR. Seems too experimental comparing with other parts of the library.

In short: I procrastinated about this. I prefer doing dumb optimizations that doesn't exercise my brain that much :P

[thehans]

We want to avoid exact convex decomposition not because the algorithm itself is hard, but because exact convex decomposition will create many parts for very simple geometries, e.g. cube - sphere.

What? cube and sphere are already convex. The convex decomposition of a convex geometry is just the one geometry

edit: ah, sorry i misinterpreted the dash as a random separator, not an difference operation.

[pca006132]

I mean cube - sphere, e.g.

[pca006132]

The algorithm I was implementing is from the paper "Offset Triangular Mesh Using the Multiple Normal Vectors of a Vertex" and "A 3D surface offset method for STL-format models". Instead of vertex normals, they use normals from incident faces/edges. Each convex edge contribute two normals, one from each incident face. For concave edges, the algorithm computes an "average" vector. Below is a diagram for the 2D case:

This is the case when there is only 1 concave edge. Note that the "normal" in this case is no longer a unit vector. When there are two consecutive concave edges, where consecutive means consider all the incident edges in cw/ccw order, there is a 3D case. However, when you get to three or more consecutive concave edges, this method doesn't really work, the paper just say something like compute the solution for every 3 consecutive normals and then average them...

After this we still have another problem: How should we compute the "blend surface" by subdivision? Again, for the convex case it is simple, just add the two vectors and normalize them. However, when we have concave edges, this does not work because the "average normal" is not a unit vector, normalizing the sum will lose the length value. Take the average of the length of the two vectors is not a solution either. I experimentally find out (read: guess and try to compute it) that if the "average normal" is the "average" of 1 or 2 normals, say $n_1 = avg(u, v), n_2 = avg(u', v')$ (where avg is the mysterious average operator, whatever that means), we can do $n_3 = avg(normalize(u, u'), normalize(v, v'))$ and $n_3$ is the vertex in the blended surface. I checked this to work with some pyramid with star-shaped base. But this doesn't work for 3 or more normals!

It would be great if someone can give me some help on the math things. I don't know much about the math here.

[zalo]

I think there's a simpler (if slower) solution to offsetting that preserves manifoldness and adds fillets...

def offset(a, distance, resolution=20):
    a_mesh = a.to_mesh()
    sphere = manifold3d.Manifold.sphere(abs(distance), resolution)
    # Hull B for each triangle of A
    a_combo = a # Should be a copy
    for tri in a_mesh.tri_verts:
        hull = manifold3d.Manifold.batch_hull(
            [sphere.translate(a_mesh.vert_properties[tri[0], :3]),
             sphere.translate(a_mesh.vert_properties[tri[1], :3]),
             sphere.translate(a_mesh.vert_properties[tri[2], :3])])

        # Union or Difference depending on the Offset Amount
        if distance > 0:
            a_combo += hull
        else:
            a_combo -= hull
    return a_combo

Original Shape:

Offset Positive took 1418.94 ms:

Offset negative took 1488.11 ms:

Minkowski Sums are a generalization of this:

def minkowski(a, b):
    a_hull    = a.hull(); b_hull = b.hull()
    a_convex  = a.get_volume() == a_hull.get_volume()
    b_convex  = b.get_volume() == b_hull.get_volume()
    a_mesh = a.to_mesh()

    # Easy Case: Just hull b for each vertex of a
    if a_convex and b_convex:
        return manifold3d.Manifold.batch_hull(
            [b.translate(v) for v in a_mesh.vert_properties[:, :3]])
    
    # Harder Case: Hull B for each triangle of A
    # Do some kind of convex decomposition here?
    if (not a_convex) and b_convex:
        a_union = a
        for tri in a_mesh.tri_verts:
            hull = manifold3d.Manifold.batch_hull(
                [b.translate(a_mesh.vert_properties[tri[0], :3]),
                 b.translate(a_mesh.vert_properties[tri[1], :3]),
                 b.translate(a_mesh.vert_properties[tri[2], :3])])
            a_union += hull
        return a_union

    # Else, we need to do a more complicated Minkowski sum with Convex Decomposition
    # Or a really really slow TriangleXTriangle Super Sum...
    assert False, "Non-Convex on Non-Convex minkowski Not Implemented Yet! Fall back to CGAL or something..."

Instead of doing it per-triangle, you’d want to do it per convex sub-component, which is where a fast exact convex decomposition would come in handy…

[elalish]

I have to admit, that brute force approach to offsetting has some nice robustness properties. I'm surprised it's that fast - obviously not great at scale, but I think most people using an offset probably want to apply it to fairly simple objects anyway. I don't think we need to do a general Minkowski op - we got started here because we noticed that the only thing anyone ever used a Minkowski for was a sphere anyway, so better to stick with offsetting for simplicity. It also gives @pca006132's CSG tree optimizations a chance to shine, with all those repeated ops.

[zalo]

It might be nice to do Linear Sweeps too. I wrote a small sweeping implementation using CoACD to see how bad it would be:

# Linear Sweep
import coacd
a_mesh = a.to_mesh()
parts = coacd.run_coacd(coacd.Mesh(a_mesh.vert_properties[:, :3], a_mesh.tri_verts)) # a list of convex hulls.
print("COACD found", len(parts), "parts")
sweep = manifold3d.Manifold()
for i, part in enumerate(parts):
    manifold_part = manifold3d.Manifold.from_mesh(manifold3d.Mesh(part[0], part[1]))
    swept_convex_part = manifold3d.Manifold.batch_hull(
        [manifold_part,
         manifold_part.translate(0.4, 0.4, 0.4)])
    sweep += swept_convex_part
return sweep

Problem here is:

[2023-11-22 10:44:45.680] [CoACD] [info] Compute Time: 4.466729499981739s
[2023-11-22 10:44:45.681] [CoACD] [info] # Convex Hulls: 72
Linear Sweep took 4584.762300015427 ms

CoACD takes 97.5% of the time 😅

Sweeping each triangle individually:

a_mesh = a.to_mesh()
swept = manifold3d.Manifold()
for tri in a_mesh.tri_verts:
    hull = manifold3d.Manifold.hull_points(
        [a_mesh.vert_properties[tri[0], :3],
         a_mesh.vert_properties[tri[1], :3],
         a_mesh.vert_properties[tri[2], :3],
         a_mesh.vert_properties[tri[0], :3] + 0.4,
         a_mesh.vert_properties[tri[1], :3] + 0.4,
         a_mesh.vert_properties[tri[2], :3] + 0.4])
    swept += hull
return swept

Results in a better mesh, and completes in only 382.6230999547988 ms

Manifold is too much faster than the supporting libraries 😂

[zalo]

@pca006132

After this we still have another problem: How should we compute the "blend surface" by subdivision? Again, for the convex case it is simple, just add the two vectors and normalize them. However, when we have concave edges, this does not work because the "average normal" is not a unit vector, normalizing the sum will lose the length value. Take the average of the length of the two vectors is not a solution either. I experimentally find out (read: guess and try to compute it) that if the "average normal" is the "average" of 1 or 2 normals, say n1=avg(u,v),n2=avg(u′,v′) (where avg is the mysterious average operator, whatever that means), we can do n3=avg(normalize(u,u′),normalize(v,v′)) and n3 is the vertex in the blended surface. I checked this to work with some pyramid with star-shaped base. But this doesn't work for 3 or more normals!

It would be great if someone can give me some help on the math things. I don't know much about the math here.

CAD Kernels (like OpenCASCADE) mesh all surfaces by subdivision of some 2D function parameterized by u,v.

For a 3-way fillet blend, I think CAD kernels try to use circular/spherical blends when possible, but will also use a polynomial to construct an interpolating surface when circular blends are not possible… I’d have to check how OpenCascade does it…

My naive recommendation would be to use barycentric interpolation of the triangle positions/normals to perform the interpolation… but I’m on the go right now, so I can’t review my usual sources 😅

If you want to see how a CAD Kernel handles it really quick, you can use my little CAD playground webapp to create a box, subtract a box from that, and select the inner three corners with the fillet tool… https://leapmotion.github.io/LeapShape/ (hopefully it all still works…)

[zalo]

One quality boost that just occurred to me is to rotate the hull spheres so that one of their vertices align to the mesh’s surface normal.

This should be a free change that addresses the flat inner face seen here without increasing the resolution:

Perhaps it can be an optional arg on the offset() function 😄

I would still like to figure out a reasonable full minkowski sum implementation though…

[elalish]

Interesting idea. Out of curiosity, what makes you attracted to the full minkowski sum op? I was always surprised it was included in OpenSCAD at all, since in all my time using that program I could never find a useful mental model for that op, so I pretty much never used it.

[zalo]

I’ll admit, the only non-convex/non-convex sum I’ve ever wanted to do is a toroid sweep along a spline path (and there are far more efficient ways of doing that).

My primary fascination with Minkowski sums has been in collision detection, where two arbitrary objects collide when the “Minkowski Difference” intersects the origin:

The bottleneck for making physics engines that don’t rely on convex primitives has always been the difficulty in computing the Minkowski sums in a reasonably fast manner (which, of course, we’ll never do exactly on meshes quickly). I find the visualization mostly educational 😄

Maybe it’s fine to have a brute force implementation for it, since few will ever use it outside of as a curiosity…

Either way, I think the convex/convex case, and the non-convex/convex case have enough fun uses that optimized implementations would make a worthwhile addition to the toolkit.

The unionRound implementation is a fun example of that:

https://github.com/TLC123/OpenSCAD_stuff

(This function in particular would benefit the most from @pca006132 ’s fast mesh offsetting method, since the slowest part is constructing a bunch of Minkowski() sphere-dilated offset meshes 😄)

[pca006132]

Interesting. I actually asked my friends doing research in robotics about usage of minkowski sum, they said that it can be used in collision detection but it is so slow that nobody really uses it.

Regarding CoCAD: It does approximate convex decomposition. I think it is faster than exact convex decomposition and give you a lot less parts. However, those approximate convex parts are not suitable for doing minkowski sum, you will get self-intersections. And it is surprising to know that the decomposition step is slower than doing hundreds of unions. I guess our optimizations are worth it.

Btw thanks for the pointer for computing blends, will look into it. I am lacking a lot of the basic knowledge about graphics and CAD.

[pca006132]

I tested the per triangle minkowski sum with openscad, it is much slower than convex decomposition when the operand is large, e.g. spheres with $fn=100, which a lot of users do use.

[zalo]

I think approximate methods are slower because they are iterative, and aiming for the smallest number of hulls possible.

I’m hopeful there are fast and exact methods (like this guys’s greedy one: https://www.reddit.com/r/GraphicsProgramming/s/0jGSbAPnxF ) that still produce vastly fewer convex parts than the number of triangles on the mesh 😄

I’d probably try to implement it by modifying the existing hull algorithm…

Could you share the OpenSCAD code you tested and the relative timings?

[zalo]

I have a new dumb algorithm in mind for convex decomposition:

• Find the most concave edge in the manifold
• Split the manifold by a plane along this edge.
• Decompose the resulting manifolds and repeat until there are no more concave edges.

It could be interesting to select the plane’s splitting angle, there are a few options:

• Pick the bisecting plane
• Pick a plane parallel to one of the adjoining triangles
• Pick a plane parallel to both adjoining triangles
• Pick a plane snapped to the nearest cardinal direction (may simplify the resulting chunks?)

My cheater heuristic will probably be:

• If the concave angle is 90°, then cut parallel to each of the adjoining triangles
• Else, cut along the bisector. Will need to test this on a bunch of meshes (both organic and inorganic), but I think this might give the best cuts on the most meshes.

And, depending on how fast “split by plane” is, this could all run at usable speeds 😄

The CGAL Pages describing their algorithms have been wonderful for clarifying my thinking: https://doc.cgal.org/latest/Convex_decomposition_3/index.html https://doc.cgal.org/latest/Minkowski_sum_3/index.html

Looking forward to trying this out once I’m back from holiday 😄

[pca006132]

I feel like picking a plane snapped to the nearest cardinal direction is probably the fastest way of split by plane, considering how the collider works.

I was thinking about something similar for the algorithm I was working on, but instead of finding concave edges I try to find pairs of faces that the swept volume will intersect (which I call conflicting faces). With those pairs, just do the decomposition by the bisecting plane. For efficiency, try to randomize a bit and try several possible cuts for each level (a bit like Monte Carlo tree search but not as sophisticated), as one cut may eliminate multiple pairs of conflicting faces. It will also prioritize cuts that have roughly equal number of pairs of conflicting faces on both parts, as we can do the decomposition of two parts in parallel.

This can work for the algorithm I was working on because it does not require the parts to be convex, this is just to avoid having conflicting faces.

And yeah the CGAL algorithm description is amazing.

[elalish]

Interesting. These ideas remind me a bit of binary space partitioning trees (BSP), which are sometimes used for mesh booleans. Sounds like they might be better for this.

[zalo]

Bisecting plane sounds good, but how do you calculate it for the swept surface?

Interesting. These ideas remind me a bit of binary space partitioning trees (BSP), which are sometimes used for mesh booleans. Sounds like they might be better for this.

(Author of RealTimeCSG and Chisel, the best brush package for Unity)

I think the tricky part is that going to brushes from meshes needs a convex decomposition 😄

I’m tempted to add a “method” arg to the function for the different speed and exactness tradeoffs… it would be cool to be able to trade exactness for not having your computer crash trying to do the operation…

I think a good additional heuristic for splitting planes is trying to get it to contain multiple concave edges at once, which should vastly cut down on the number of convex parts for “blocky” objects (though it’ll do little for “organic” objects).

[pca006132]

No, I am I only calculate it for two planes, because I only need to separate them so their swept volumes will not be intersecting.

For the heuristic, I think apart from whatever geometry things that you come up with, adding Monte Carlo tree search (e.g. try not only the edge with the highest weight, but several ones and use the best one) will probably help a bit.

[zalo]

Hrmmm, testing out my Solid Exact Convex Decomposition, it seems like it takes about as long as CoACD and VHACD, but outputs 2-6x more convex parts (close to the order of as many concave edges as there are on the model)...

So far, not seeing the benefit except for very blocky models, which it decomposes great.

Perhaps it's time to explore the surface decomposition methods... https://dpd.cs.princeton.edu/Papers/ChazelleDobkinShourabouraTal.pdf (Which I'm skeptical of, since the pieces could generate hulls that extend outside the bounds of the original part...)

Or the CGAL Method: https://alexandria.tue.nl/openaccess/Metis212139.pdf

Or perhaps another one...

[pca006132]

The issue with the number of convex parts is the reason we think that mesh offsetting using minkowski sum and convex decomposition will not scale. Doesn't seem solvable without doing some approximation, but existing approximate convex decompositions are not suitable for minkowski sum.

For the openscad test, I will post it later this week. I am currently in hospital and does not have access to my computer for several days. I was just modifying the minkowski sum code for manifold backend in openscad, and run them against some benchmarks (in the manifold integration PR for openscad).

[pca006132]

Btw, wondering if porting CGAL's convex decomposition function will make it significantly faster. For openscad, iirc the bottleneck is currently in the convex decomposition part.

CGAL is only accepts Nef polyhedron for the convex decomposition, which uses exact arithmetic and operations on them are known to be pretty slow. They don't have any parallelization for the decomposition, hopefully we can find some opportunities for that.

[elalish]

Regarding @zalo's simple offset method - does negative offsetting work when you use convex decomposition? That feels like a pretty useful feature, but I'm not sure that Minkowski can do that. I would prefer a spherical offset that goes negative and positive, to a Minkowski that takes different convex objects but only goes one way.

[pca006132]

For negative offset, just subtract the object from the bounding box, do a positive offset, and subtract again?