|
| 1 | +import math |
| 2 | +import compas_libigl as igl |
| 3 | +from compas.colors import Color, ColorMap |
| 4 | +from compas.geometry import Plane |
| 5 | +from compas.datastructures import Mesh |
| 6 | +from compas.geometry import Rotation |
| 7 | +from compas.geometry import Scale |
| 8 | +from compas_viewer import Viewer |
| 9 | + |
| 10 | + |
| 11 | +def compute_plane_distance(point, plane): |
| 12 | + """Compute signed distance from a point to a plane.""" |
| 13 | + vector = plane.point - point |
| 14 | + return plane.normal.dot(vector) |
| 15 | + |
| 16 | + |
| 17 | +def slice_mesh_with_planes(mesh, planes, attribute_name="cut", combine="min"): |
| 18 | + """Slice a mesh using the combined distance field from multiple planes.""" |
| 19 | + if not planes: |
| 20 | + return None, None |
| 21 | + |
| 22 | + # Compute signed distances for all planes |
| 23 | + all_distances = [] |
| 24 | + for plane in planes: |
| 25 | + distances = [] |
| 26 | + for vertex in mesh.vertices(): |
| 27 | + point = mesh.vertex_attributes(vertex, "xyz") |
| 28 | + dist = compute_plane_distance(point, plane) |
| 29 | + distances.append(dist) |
| 30 | + all_distances.append(distances) |
| 31 | + |
| 32 | + # Combine distances based on the chosen method |
| 33 | + if len(all_distances) == 1: |
| 34 | + combined_distances = all_distances[0] |
| 35 | + else: |
| 36 | + combined_distances = all_distances[0] # Start with first plane's distances |
| 37 | + for distances in all_distances[1:]: # Combine with remaining planes |
| 38 | + if combine == "min": |
| 39 | + combined_distances = [min(d1, d2) for d1, d2 in zip(combined_distances, distances)] |
| 40 | + elif combine == "max": |
| 41 | + combined_distances = [max(d1, d2) for d1, d2 in zip(combined_distances, distances)] |
| 42 | + elif combine == "average": |
| 43 | + combined_distances = [(d1 + d2) / 2 for d1, d2 in zip(combined_distances, distances)] |
| 44 | + elif combine == "multiply": |
| 45 | + combined_distances = [d1 * d2 for d1, d2 in zip(combined_distances, distances)] |
| 46 | + |
| 47 | + for vertex, distance in zip(mesh.vertices(), combined_distances): |
| 48 | + mesh.vertex_attribute(vertex, attribute_name, distance) |
| 49 | + |
| 50 | + # Remesh along the zero isoline (intersection) |
| 51 | + V2, F2, L = igl.trimesh_remesh_along_isoline( |
| 52 | + mesh.to_vertices_and_faces(), |
| 53 | + mesh.vertices_attribute(attribute_name), |
| 54 | + 0 |
| 55 | + ) |
| 56 | + |
| 57 | + # Split faces based on labels |
| 58 | + below_faces = [f for i, f in enumerate(F2) if L[i] == 0] |
| 59 | + above_faces = [f for i, f in enumerate(F2) if L[i] == 1] |
| 60 | + |
| 61 | + # Get unique vertices for each part |
| 62 | + below_vertices = set() |
| 63 | + above_vertices = set() |
| 64 | + for face in below_faces: |
| 65 | + below_vertices.update(face) |
| 66 | + for face in above_faces: |
| 67 | + above_vertices.update(face) |
| 68 | + |
| 69 | + # Create vertex maps for new indices |
| 70 | + below_vmap = {old: new for new, old in enumerate(sorted(below_vertices))} |
| 71 | + above_vmap = {old: new for new, old in enumerate(sorted(above_vertices))} |
| 72 | + |
| 73 | + # Create new vertex lists |
| 74 | + below_verts = [V2[i] for i in sorted(below_vertices)] |
| 75 | + above_verts = [V2[i] for i in sorted(above_vertices)] |
| 76 | + |
| 77 | + # Remap face indices |
| 78 | + below_faces = [[below_vmap[v] for v in face] for face in below_faces] |
| 79 | + above_faces = [[above_vmap[v] for v in face] for face in above_faces] |
| 80 | + |
| 81 | + # Create new meshes |
| 82 | + below_mesh = Mesh.from_vertices_and_faces(below_verts, below_faces) if below_faces else None |
| 83 | + above_mesh = Mesh.from_vertices_and_faces(above_verts, above_faces) if above_faces else None |
| 84 | + |
| 85 | + return below_mesh, above_mesh |
| 86 | + |
| 87 | + |
| 88 | +def recursive_slice(mesh, plane_groups, depth=0, region_id=0): |
| 89 | + """Recursively slice mesh with groups of planes and assign region IDs.""" |
| 90 | + if depth >= len(plane_groups): |
| 91 | + return [(mesh, region_id)] |
| 92 | + |
| 93 | + below_mesh, above_mesh = slice_mesh_with_planes(mesh, plane_groups[depth], combine="min") |
| 94 | + |
| 95 | + results = [] |
| 96 | + # Process below mesh (region_id stays the same) |
| 97 | + if below_mesh and below_mesh.number_of_faces() > 0: |
| 98 | + results.extend(recursive_slice(below_mesh, plane_groups, depth + 1, region_id)) |
| 99 | + |
| 100 | + # Process above mesh (region_id gets modified) |
| 101 | + if above_mesh and above_mesh.number_of_faces() > 0: |
| 102 | + results.extend(recursive_slice(above_mesh, plane_groups, depth + 1, region_id + 2**depth)) |
| 103 | + |
| 104 | + return results |
| 105 | + |
| 106 | + |
| 107 | +# Load and transform mesh |
| 108 | +mesh = Mesh.from_off(igl.get_beetle()) |
| 109 | +Rx = Rotation.from_axis_and_angle([1, 0, 0], math.radians(90)) |
| 110 | +Rz = Rotation.from_axis_and_angle([0, 0, 1], math.radians(90)) |
| 111 | +S = Scale.from_factors([10, 10, 10]) |
| 112 | +mesh.transform(S * Rz * Rx) |
| 113 | + |
| 114 | +# Define different directions for planes |
| 115 | +directions = [ |
| 116 | + [0, 1, 1], # Direction 1 |
| 117 | + # [1, 0, 1], # Direction 2 |
| 118 | + # # Add more directions if needed |
| 119 | + # [1, 1, 0], # Direction 3 |
| 120 | +] |
| 121 | + |
| 122 | +# Create groups of planes at each offset |
| 123 | +plane_groups = [] |
| 124 | +for i in range(-2, 3): # Creates 5 offset positions |
| 125 | + # At each offset, create a group with planes in different directions |
| 126 | + planes_at_offset = [Plane([0, 0, i], direction) for direction in directions] |
| 127 | + plane_groups.append(planes_at_offset) |
| 128 | + |
| 129 | +# Perform recursive slicing |
| 130 | +mesh_regions = recursive_slice(mesh, plane_groups) |
| 131 | +print(f"Number of regions: {len(mesh_regions)}") |
| 132 | + |
| 133 | +# Setup visualization |
| 134 | +viewer = Viewer() |
| 135 | + |
| 136 | +# Create color gradient |
| 137 | +cmap = ColorMap.from_mpl("viridis") |
| 138 | + |
| 139 | +# Add each region with its color |
| 140 | +for mesh_part, region_id in mesh_regions: |
| 141 | + if mesh_part and mesh_part.number_of_faces() > 0: |
| 142 | + # Normalize region_id to [0,1] for coloring |
| 143 | + color = cmap(region_id / (2**len(plane_groups))) |
| 144 | + viewer.scene.add(mesh_part, facecolor=color, show_lines=False) |
| 145 | + |
| 146 | +viewer.show() |
0 commit comments