New Convex Decomposition Technique using Weighted-Voronoi Diagrams aka Power Diagrams aka Laguerre Tessellations

voronoi_decomposition.py

import tess
import random
import trimesh
import manifold3d
import numpy as np
from time import perf_counter

rand_color = [random.random(), random.random(), random.random()]
def explode(convex_pieces, explode_amount = 1.5, debug_shapes = None):
    global rand_color
    exploded_pieces = []
    for i, convex_piece in enumerate(convex_pieces):
        centroid = np.mean(convex_piece.to_mesh().vert_properties[:, :3], axis=0)
        offset = centroid*explode_amount - centroid
        exploded_piece = convex_piece.translate(offset[0], offset[1], offset[2])
        rand_color = [random.random(), random.random(), random.random()]
        try:
            exploded_piece  = exploded_piece .set_properties(3, lambda pos, oldProp: rand_color)
            if debug_shapes is not None:
                debug_shapes[i] = debug_shapes[i].set_properties(3, lambda pos, oldProp: rand_color)
                exploded_pieces.append(debug_shapes[i])
        except:
            pass
            #print("Failed to set properties")
        exploded_pieces.append(exploded_piece)
    return exploded_pieces

def save_mesh(model, name, ext="glb"):
    mesh = model.to_mesh()
    if mesh.vert_properties.shape[1] > 3:
        vertices =  mesh.vert_properties[:, :3]
        colors   = (mesh.vert_properties[:, 3:] * 255).astype(np.uint8)
    else:
        vertices = mesh.vert_properties
        colors   = None
    meshOut = trimesh.Trimesh(vertices=vertices, faces=mesh.tri_verts,
                              vertex_colors=colors)
    trimesh.exchange.export.export_mesh(meshOut, name+'.'+ext, ext)
    print('Exported model to ', name+'.'+ext)

def circumcircle(triangle_vertices):
    a            = triangle_vertices[0] - triangle_vertices[2]
    b            = triangle_vertices[1] - triangle_vertices[2]
    a_length     = np.linalg.norm(a); b_length = np.linalg.norm(b)
    numerator    = np.cross((((a_length*a_length)*b) - ((b_length*b_length)*a)), np.cross(a, b))
    crs          = np.linalg.norm(np.cross(a, b))
    denominator  = 2.0 * (crs * crs)
    circumcenter = (numerator/denominator) + triangle_vertices[2]
    circumradius = np.linalg.norm(circumcenter - triangle_vertices[0])
    return circumcenter, circumradius

debug_shapes = []
original_to_joggle = {}; joggle_to_original = {}
def convex_decomposition(shape : manifold3d.Manifold, enforce_hulls : bool = False) -> list[manifold3d.Manifold]:
    global debug_shapes, original_to_joggle, joggle_to_original
    outputs = []
    if shape is None:
        print("[ERROR] SHAPE IS NONE!!!")
        return []
    shapes = shape.decompose()
    if len(shapes) == 0:
        print("[ERROR] INVALID DECOMPOSITION!!!")
        return [shape]
    for shape in shapes:
        if shape is None:
            continue
        mesh        = shape.to_mesh()
        vertices    = (mesh.vert_properties[:, :3]).astype(np.float64)
        cur_trimesh = trimesh.Trimesh(vertices=vertices, faces=mesh.tri_verts)

        t0_python_heavy = perf_counter()
        # Get the adjacent faces of each concave segment
        nonconvex_face_pairs     = ~(cur_trimesh.face_adjacency_projections < 1e-6)
        nonconvex_adjacent_faces = cur_trimesh .face_adjacency[nonconvex_face_pairs]
        unique_faces             = np.unique(nonconvex_adjacent_faces.flatten())

        # Compute First Pass Voronoi Cells
        voronoi_cells = np.zeros((len(unique_faces), 3), dtype=np.float64)
        voronoi_radii = np.zeros( len(unique_faces),     dtype=np.float64)
        for i, face in enumerate(unique_faces):
            face_verts = vertices[cur_trimesh.faces[face]]
            voronoi_cells[i], voronoi_radii[i] = circumcircle(face_verts)

        # Wrong: Take only the unique voronoi cells
        #voronoi_cells, unique_indices = np.unique(voronoi_cells.round(decimals=3), axis=0, return_index=True)
        #voronoi_radii = voronoi_radii[unique_indices]

        # Joggle only triangles with identical circumcenters!
        _, unique_indices = np.unique(voronoi_cells.round(decimals=3), axis=0, return_index=True)
        duplicate_circ_indices = np.setdiff1d(np.arange(len(unique_faces)), unique_indices)
        original_to_joggle = {}; joggle_to_original = {}
        for i, face in enumerate(duplicate_circ_indices):
            orig_face_verts    = vertices[cur_trimesh.faces[face]].copy()[0]
            joggle_face_verts  = orig_face_verts + np.random.random(3) * 0.0000001
            original_to_joggle[  orig_face_verts.round(decimals=3).tobytes()] = (joggle_face_verts[0], joggle_face_verts[1], joggle_face_verts[2])
            joggle_to_original[joggle_face_verts.round(decimals=3).tobytes()] =   orig_face_verts
        def joggle(pos):
            global original_to_joggle
            key = np.array(pos, dtype=np.float64).round(decimals=3).tobytes()
            if key in original_to_joggle:
                #print("JOGGLING TO", original_to_joggle[key], "FROM", pos)
                return original_to_joggle[key]
            return pos
        joggled_shape    = shape.warp(joggle)
        joggled_vertices = (joggled_shape.to_mesh().vert_properties[:, :3]).astype(np.float64)
        for i, face in enumerate(unique_faces):
            face_verts = joggled_vertices[cur_trimesh.faces[face]]
            voronoi_cells[i], voronoi_radii[i] = circumcircle(face_verts)

        #for i, cell in enumerate(voronoi_cells):
        #    # DEBUG SPHERES Show the voronoi cells
        #    debug_shapes += [manifold3d.Manifold.sphere(voronoi_radii[i] * 0.1, 30).translate(cell)]

        print("Python Heavy took", (perf_counter()-t0_python_heavy)*1000, "ms")

        #print("Voronoi Cells:", voronoi_cells.shape, "Voronoi Radii:", voronoi_radii.shape, "Voronoi Cells:", voronoi_cells.dtype, "Voronoi Radii:", voronoi_radii.dtype)
        #print("Voronoi Cells:", voronoi_cells, "Voronoi Radii:", voronoi_radii)

        t0_voronoi = perf_counter()
        voronoi_diagram = tess.Container(voronoi_cells, limits=((-2,-2,-2), (2,2,2)), periodic=False, radii=voronoi_radii)
        print("Voronoi Diagram took", (perf_counter()-t0_voronoi)*1000, "ms")

        for region in voronoi_diagram:
            region_vertices = np.array(region.vertices(), dtype=np.float64)

            # Ignore the empty and infinite regions
            if len(region_vertices) == 0:
                continue

            # Regions must be at least tetrahedra
            if len(region_vertices) > 3:
                # Unjoggle the voronoi region vertices!!
                for i, vertex in enumerate(region_vertices):
                    key = vertex.round(decimals=3).tobytes()
                    if key in joggle_to_original:
                        region_vertices[i] = joggle_to_original[key]

                cur_hull  = manifold3d.Manifold.hull_points(region_vertices)

                cur_hull ^= shape

                if cur_hull.get_volume() > 0.0:
                    outputs  += [cur_hull] if not enforce_hulls else [manifold3d.Manifold.hull(cur_hull)]
    return outputs

if __name__ == "__main__":
    sphere = manifold3d.Manifold.sphere(0.6, 20)
    cube   = manifold3d.Manifold.cube  (1.0, 1.0, 1.0, True)

    #fun_shape = cube - sphere
    #fun_shape = cube + cube.translate(0.5, 0.5, 0.5)
    fun_shape = sphere + sphere.translate(0.5, 0.5, 0.5)
    save_mesh(fun_shape, "fun_shape")

    t0 = perf_counter()
    convex_parts      = convex_decomposition(fun_shape, True)
    print("Convex Decomposition took", (perf_counter()-t0)*1000, "ms creating", len(convex_parts), "pieces")

    # Calculate volume of original mesh and union of convex hulls
    original_volume = fun_shape.get_volume()
    convex_volume   = 0.0
    union           = manifold3d.Manifold()
    for convex_part in convex_parts:
        union += convex_part
    union_volume = union.get_volume()
    print("Original Volume:", original_volume, "Convex Volume:", union_volume, "Difference:", str(round(((union_volume - original_volume)/original_volume) * 100.0, 5))+"%")

    print("Convex Parts:", len(convex_parts), "Debug Shapes:", len(debug_shapes))
    convex_decomp = manifold3d.Manifold.compose(explode(convex_parts, 1.0005, debug_shapes))
    save_mesh(convex_decomp, "convex_decomposition", "glb")

    exploded_convex_decomp = manifold3d.Manifold.compose(explode(convex_parts, 1.5))
    save_mesh(exploded_convex_decomp, "exploded_convex_decomposition", "glb")

    exploded_convex_decomp = manifold3d.Manifold.compose(explode(convex_parts, 1.5))
    save_mesh(exploded_convex_decomp, "exploded_convex_decomposition", "stl")

    #voronoi_debug = manifold3d.Manifold.compose([fun_shape] + debug_shapes)
    #save_mesh(voronoi_debug, "voronoi_debug", "stl")