offset a convex manifold

[yxdragon]

I am working on discrete element method (DEM) particle generation. I need to offset(dilation or erosion) a manifold. as I saw your discussion about how to offset a arbitrary manifold is hard. Now my manifolds are all convex. So what is the best way to inflate or dilate?

1. I can make a huge enough cubic, then count every face's normal vector, then use trim_by_plane.
2. about the inflation in round mode, I build a sphere, and extract all verts, then I use numpy to broadcast the sphere on the manifold's every vert. then make a convexhull.

It seems that, count vert's normal vector then offset the verts is more efficient than trim by every plane. and it is ok for dilation. but when erosion, some face may be "swallowed", and this method not works.

Is it ok? or is there any better way?

[yxdragon]

I think manifold3d could provide a method, generate a manifold by the given planes (normal vector + distance), just like convex hull. I found it is a linear program problem.

[pca006132]

I think for erosion, the typical way is to subtract the mesh from its bounding box, dilate the subtraction result, and perform the subtraction again.

For generating manifold by the given planes, maybe you can try SDF?

[elalish]

If you look through #666 you'll find a bunch of related work and discussions, along with links to some other attempted PRs. It's a hard problem in general. For a convex object, the convex hull is probably the way to go. Vert offset will not be accurate even when it doesn't swallow faces.

[yxdragon]

I think for erosion, the typical way is to subtract the mesh from its bounding box, dilate the subtraction result, and perform the subtraction again.

For generating manifold by the given planes, maybe you can try SDF?

Yes, I implemented it by "cut face by face" as you said. but I think the trim_by_plane's perfomance would slow down, as the bounding cube's face growing.

And I had tested it is true. So I try another way: trim the initial cube by every face plan, then run a batch_bool. It is 3 times faster.

def convex_offset(mani, distance):
    mesh = mani.to_mesh()
    x1, y1, z1, x2, y2, z2 = mani.bounding_box()
    d2 = abs(distance) * 3
    box = mfd.Manifold.cube([(x2-x1)+d2,(y2-y1)+d2,(z2-z1)+d2])
    box = box.translate([x1-d2/2, y1-d2/2, z1-d2/2])
    verts, faces = mesh.vert_properties, mesh.tri_verts
    p0, p1, p2 = verts[faces[:,0]], verts[faces[:,1]], verts[faces[:,2]]
    v1, v2 = p1 - p0, p1 - p2
    nv = np.cross(v1, v2);
    nv /= np.linalg.norm(nv, axis=1)[:,None]
    ds = (p1 * nv).sum(axis=1)
    rst = [box.trim_by_plane(n, s - distance) for n,s in zip(nv, ds)]
    return mfd.Manifold.batch_boolean(rst, mfd.OpType.Intersect)

[yxdragon]

import manifold3d as mfd
from time import time
import numpy as np

def mesh_ravel(verts, faces):
    verts = np.concatenate((verts, verts[:1]+np.nan), axis=0)
    faces = np.concatenate((faces, faces[:,:2]), axis=1)
    faces[:,-1] = -1; faces = faces.ravel()[:-1]
    return verts[faces]

def show_mesh(mesh, color=(1,1,1)):
    verts, faces = mesh.vert_properties, mesh.tri_verts
    polygons = mesh_ravel(verts, faces)
    from vispy.scene import SceneCanvas
    from vispy.scene.visuals import Line
    from vispy import app
    
    dim = polygons.shape[1]
    if dim==2:
        polygons = np.concatenate((polygons, polygons[:,:1]*0), axis=1)
    canvas = SceneCanvas(keys='interactive', title='Polygon', show=True)
    v = canvas.central_widget.add_view()
    v.bgcolor = (0,0,0)
    v.camera = 'panzoom' if dim==2 else 'turntable'
    v.camera.aspect = 1
    x1, y1, z1 = np.nanmin(polygons, axis=0)
    x2, y2, z2 = np.nanmax(polygons, axis=0)
    if dim==2: v.camera.zoom(((x2-x1)**2+(y2-y1)**2)**0.5)
    v.camera.center = ((x1+x2)/2, (y1+y2)/2, (z1+z2)/2)
    Line(polygons, color=color, parent=v.scene)
    app.run()

def convex_offset(mani, distance):
    mesh = mani.to_mesh()
    x1, y1, z1, x2, y2, z2 = mani.bounding_box()
    d2 = abs(distance) * 3
    box = mfd.Manifold.cube([(x2-x1)+d2,(y2-y1)+d2,(z2-z1)+d2])
    box = box.translate([x1-d2/2, y1-d2/2, z1-d2/2])
    verts, faces = mesh.vert_properties, mesh.tri_verts
    p0, p1, p2 = verts[faces[:,0]], verts[faces[:,1]], verts[faces[:,2]]
    v1, v2 = p1 - p0, p1 - p2
    nv = np.cross(v1, v2);
    nv /= np.linalg.norm(nv, axis=1)[:,None]
    ds = (p1 * nv).sum(axis=1)
    rst = [box.trim_by_plane(n, s - distance) for n,s in zip(nv, ds)]
    return mfd.Manifold.batch_boolean(rst, mfd.OpType.Intersect)

def dilateion_round(mani, distance, accu=8):
    sphere = mfd.Manifold.sphere(distance, accu)
    pts = sphere.to_mesh().vert_properties
    verts = mani.to_mesh().vert_properties
    verts = (verts[:,None,:] + pts).reshape(-1,3)
    return mfd.Manifold.hull_points(verts)

def offset3d_cvx(mani, distance, jointype=mfd.JoinType.Round, accu=8):
    if distance < 0 or jointype != mfd.JoinType.Round:
        return convex_offset(mani, distance)
    return dilateion_round(mani, distance, accu)


if __name__ == '__main__':
    cube = mfd.Manifold.cube([3,3,3])

    dilate = offset3d_cvx(cube, 0.2)
    erode = offset3d_cvx(cube, -0.5)
    both = mfd.Manifold.compose([dilate, erode])

    show_mesh(both.to_mesh(), (1,1,1))

[yxdragon]

it 's also could be used as a smooth method.
offset -r, then offset r with round mode.

[pca006132]

Yes, openscad users often use this to smooth things out, but convex decomposition is expensive.

[yxdragon]

def convex_offset(mani, distance):
    mesh = mani.to_mesh()
    x1, y1, z1, x2, y2, z2 = mani.bounding_box()
    d2 = abs(distance) * 3
    box = mfd.Manifold.cube([(x2-x1)+d2,(y2-y1)+d2,(z2-z1)+d2])
    box = box.translate([x1-d2/2, y1-d2/2, z1-d2/2])
    verts, faces = mesh.vert_properties, mesh.tri_verts
    p0, p1, p2 = verts[faces[:,0]], verts[faces[:,1]], verts[faces[:,2]]
    v1, v2 = p1 - p0, p1 - p2
    nv = np.cross(v1, v2);
    nv /= np.linalg.norm(nv, axis=1)[:,None]
    ds = (p1 * nv).sum(axis=1)
    rst = [box.trim_by_plane(n, s - distance) for n,s in zip(nv, ds)]
    return mfd.Manifold.batch_boolean(rst, mfd.OpType.Intersect)

this function is tested ok, but not fast enough. can manifold3d give a C+ implementation?

[pca006132]

C++ version may not be much faster, it is very likely that the majority of time is spent within batch boolean, and we can't do much to optimize it further for now.

Did you try broadcasting the points and do convex hull? Will it be faster?

[yxdragon]

"broadcasting points" just works for "dilation with round mode", and it has the same performance as batch boolean. (but is faster than trim a cube face by face)

[pca006132]

Anyway, I don't think there is much we can do here for now, and we don't plan to support this only for convex manifold. We should continue discussion in #192 or open some new discussions. Closing now.

[yxdragon]

@elalish does the mesh.face_id is credible?

def convex_offset(mani, distance):
    mesh = mani.to_mesh()
    x1, y1, z1, x2, y2, z2 = mani.bounding_box()
    d2 = abs(distance) * 3
    box = mfd.Manifold.cube([(x2-x1)+d2,(y2-y1)+d2,(z2-z1)+d2])
    box = box.translate([x1-d2/2, y1-d2/2, z1-d2/2])
    verts, faces = mesh.vert_properties, mesh.tri_verts
    p0, p1, p2 = verts[faces[:,0]], verts[faces[:,1]], verts[faces[:,2]]
    v1, v2 = p1 - p0, p1 - p2
    nv = np.cross(v1, v2);
    nv /= np.linalg.norm(nv, axis=1)[:,None]
    ds = (p1 * nv).sum(axis=1)
    rst = [box.trim_by_plane(n, s - distance) for n,s in zip(nv, ds)]
    return mfd.Manifold.batch_boolean(rst, mfd.OpType.Intersect)

I use trim_by_plan face by face to offset a convex manifold. but I found that I need not cut by every tri face, (and it generate many small segment in a face for cutting by coplanar tri face many times)

So I use face_id to get the first unique face to cut.

 _, idx = np.unique(mesh.face_id, return_index=True)
verts, faces = mesh.vert_properties, mesh.tri_verts[idx]

and the result is ok in most case. but with some bug.

[elalish]

If you think you've found a bug, please make a minimal repro - ideally in the form of a PR that adds a TEST that doesn't pass.