This tool is particularly useful for creating closed 3D printable meshes from one or several open meshes.
I made it using ChatGPT and it took a couple of hours of optimizing and testing.
It can be used as an alternative to ShrinkWrap.
Initial input, open mesh:
Voxel size = 5mm
Total computing time = 3.56 sec
Voxel size = 1mm
Total computing time = 41.51 sec
Voxel size = 0.5mm
Total computing time = 203.36 sec
# VOXEL TOOL V3
import Rhino
import rhinoscriptsyntax as rs
import scriptcontext as sc
import System
import math
import time
from concurrent.futures import ThreadPoolExecutor
from System.Collections.Generic import List
# ================================================================
# SETTINGS
# ================================================================
DEFAULT_VOXEL_SIZE = 100.0
WORKERS = 5
# Sampling spacing relative to voxel size.
# 0.50 = reliable
SAMPLE_FACTOR = 0.50
# Number of voxel expansion layers
DILATION_STEPS = 1
# Bounding-box padding
PADDING_VOXELS = 2
# ================================================================
# INPUT
# ================================================================
def get_mesh_ids():
ids = rs.SelectedObjects()
if ids:
result = []
for obj_id in ids:
obj = sc.doc.Objects.FindId(obj_id)
if (
obj and
isinstance(
obj.Geometry,
Rhino.Geometry.Mesh
)
):
result.append(obj_id)
if result:
return result
return rs.GetObjects(
"Select meshes to voxelize",
rs.filter.mesh,
preselect=True
)
# ================================================================
# COMBINE INPUT MESHES
# ================================================================
def combine_meshes(ids):
combined = Rhino.Geometry.Mesh()
for obj_id in ids:
obj = sc.doc.Objects.FindId(obj_id)
if not obj:
continue
mesh = obj.Geometry
if not isinstance(
mesh,
Rhino.Geometry.Mesh
):
continue
combined.Append(mesh)
combined.Vertices.CullUnused()
# Temporary triangulation only.
# Original Rhino meshes are untouched.
combined.Faces.ConvertQuadsToTriangles()
combined.Compact()
return combined
# ================================================================
# VOXEL INTEGER KEYS
# ================================================================
def voxel_key(i, j, k, nx, ny):
return i + nx * (j + ny * k)
# ================================================================
# EXTRACT TRIANGLES
# ================================================================
def extract_triangles(mesh):
vertices = mesh.Vertices
triangles = []
append = triangles.append
for index in range(mesh.Faces.Count):
f = mesh.Faces[index]
a = vertices[f.A]
b = vertices[f.B]
c = vertices[f.C]
append(
(
float(a.X),
float(a.Y),
float(a.Z),
float(b.X),
float(b.Y),
float(b.Z),
float(c.X),
float(c.Y),
float(c.Z)
)
)
return triangles
# ================================================================
# TRIANGLE RASTERIZER
# ================================================================
def rasterize_chunk(
triangles,
voxel_size,
xmin,
ymin,
zmin,
nx,
ny,
nz
):
occupied = set()
add_voxel = occupied.add
sample_spacing = (
voxel_size *
SAMPLE_FACTOR
)
inv_voxel = 1.0 / voxel_size
for tri in triangles:
ax = tri[0]
ay = tri[1]
az = tri[2]
bx = tri[3]
by = tri[4]
bz = tri[5]
cx = tri[6]
cy = tri[7]
cz = tri[8]
# --------------------------------------------------------
# EDGE VECTORS
# --------------------------------------------------------
abx = bx - ax
aby = by - ay
abz = bz - az
acx = cx - ax
acy = cy - ay
acz = cz - az
bcx = cx - bx
bcy = cy - by
bcz = cz - bz
# --------------------------------------------------------
# EDGE LENGTHS
# --------------------------------------------------------
lab = math.sqrt(
abx * abx +
aby * aby +
abz * abz
)
lac = math.sqrt(
acx * acx +
acy * acy +
acz * acz
)
lbc = math.sqrt(
bcx * bcx +
bcy * bcy +
bcz * bcz
)
longest = max(
lab,
lac,
lbc
)
divisions = max(
1,
int(
math.ceil(
longest /
sample_spacing
)
)
)
inv_div = 1.0 / divisions
# ========================================================
# BARYCENTRIC SAMPLING
# ========================================================
for ui in range(divisions + 1):
u = ui * inv_div
remaining = divisions - ui
for vi in range(remaining + 1):
v = vi * inv_div
w = 1.0 - u - v
x = (
w * ax +
u * bx +
v * cx
)
y = (
w * ay +
u * by +
v * cy
)
z = (
w * az +
u * bz +
v * cz
)
i = int(
math.floor(
(x - xmin) *
inv_voxel
)
)
j = int(
math.floor(
(y - ymin) *
inv_voxel
)
)
k = int(
math.floor(
(z - zmin) *
inv_voxel
)
)
if (
0 <= i < nx and
0 <= j < ny and
0 <= k < nz
):
add_voxel(
i +
nx *
(
j +
ny * k
)
)
return occupied
# ================================================================
# MULTITHREADED TRIANGLE RASTERIZATION
# ================================================================
def rasterize_triangles(
triangles,
voxel_size,
xmin,
ymin,
zmin,
nx,
ny,
nz
):
count = len(triangles)
worker_count = min(
WORKERS,
count
)
if worker_count < 1:
return set()
chunk_size = int(
math.ceil(
float(count) /
worker_count
)
)
chunks = []
for start in range(
0,
count,
chunk_size
):
chunks.append(
triangles[
start:
start + chunk_size
]
)
print("")
print(
"Rasterizing with {} workers...".format(
len(chunks)
)
)
occupied = set()
with ThreadPoolExecutor(
max_workers=worker_count
) as executor:
futures = []
for chunk in chunks:
futures.append(
executor.submit(
rasterize_chunk,
chunk,
voxel_size,
xmin,
ymin,
zmin,
nx,
ny,
nz
)
)
for index, future in enumerate(futures):
result = future.result()
occupied.update(result)
print(
"Worker {} / {} complete "
"surface voxels: {:,}".format(
index + 1,
len(futures),
len(occupied)
)
)
return occupied
# ================================================================
# FAST 6-NEIGHBOUR DILATION
# ================================================================
def dilate_voxels(
occupied,
nx,
ny,
nz,
steps
):
if steps <= 0:
return occupied
print("")
print(
"Expanding voxel shell by {} layer(s)...".format(
steps
)
)
current = occupied
layer = nx * ny
for step in range(steps):
expanded = set(current)
add = expanded.add
for key in current:
k = key // layer
rem = key - k * layer
j = rem // nx
i = rem - j * nx
if i > 0:
add(key - 1)
if i < nx - 1:
add(key + 1)
if j > 0:
add(key - nx)
if j < ny - 1:
add(key + nx)
if k > 0:
add(key - layer)
if k < nz - 1:
add(key + layer)
current = expanded
print(
" layer {}: {:,} voxels".format(
step + 1,
len(current)
)
)
return current
# ================================================================
# OPTIMIZED OUTPUT MESH
#
# Main improvements:
#
# - final mesh vertex indices are assigned DURING boundary scan
# - no used_vertices set
# - no second vertex-map construction
# - faces already contain final mesh vertex indices
# - Point3f instead of Point3d
# - bulk AddVertices()
# - bulk AddFaces()
# ================================================================
def build_voxel_mesh_fast(
occupied,
xmin,
ymin,
zmin,
voxel_size,
nx,
ny,
nz
):
print("")
print("Building optimized boundary mesh...")
total_start = time.time()
# Grid vertex dimensions
gvx = nx + 1
gvy = ny + 1
gv_layer = gvx * gvy
voxel_layer = nx * ny
# ------------------------------------------------------------
# GRID VERTEX KEY -> FINAL MESH INDEX
# ------------------------------------------------------------
vertex_map = {}
# Store only the grid keys in final vertex-index order.
vertex_keys = []
vertex_append = vertex_keys.append
# Faces already reference final mesh indices.
face_indices = []
face_append = face_indices.append
contains = occupied.__contains__
get_existing = vertex_map.get
# ------------------------------------------------------------
# FAST VERTEX INDEX CREATION
# ------------------------------------------------------------
def get_vertex_index(ix, iy, iz):
grid_key = (
ix +
gvx *
(
iy +
gvy * iz
)
)
existing = get_existing(grid_key)
if existing is not None:
return existing
index = len(vertex_keys)
vertex_map[grid_key] = index
vertex_append(grid_key)
return index
# ============================================================
# SINGLE BOUNDARY SCAN
# ============================================================
boundary_start = time.time()
for key in occupied:
k = key // voxel_layer
rem = key - k * voxel_layer
j = rem // nx
i = rem - j * nx
# ========================================================
# -X
# ========================================================
if (
i == 0 or
not contains(key - 1)
):
a = get_vertex_index(
i,
j,
k
)
b = get_vertex_index(
i,
j,
k + 1
)
c = get_vertex_index(
i,
j + 1,
k + 1
)
d = get_vertex_index(
i,
j + 1,
k
)
face_append(
(a, b, c, d)
)
# ========================================================
# +X
# ========================================================
if (
i == nx - 1 or
not contains(key + 1)
):
a = get_vertex_index(
i + 1,
j,
k
)
b = get_vertex_index(
i + 1,
j + 1,
k
)
c = get_vertex_index(
i + 1,
j + 1,
k + 1
)
d = get_vertex_index(
i + 1,
j,
k + 1
)
face_append(
(a, b, c, d)
)
# ========================================================
# -Y
# ========================================================
if (
j == 0 or
not contains(key - nx)
):
a = get_vertex_index(
i,
j,
k
)
b = get_vertex_index(
i + 1,
j,
k
)
c = get_vertex_index(
i + 1,
j,
k + 1
)
d = get_vertex_index(
i,
j,
k + 1
)
face_append(
(a, b, c, d)
)
# ========================================================
# +Y
# ========================================================
if (
j == ny - 1 or
not contains(key + nx)
):
a = get_vertex_index(
i,
j + 1,
k
)
b = get_vertex_index(
i,
j + 1,
k + 1
)
c = get_vertex_index(
i + 1,
j + 1,
k + 1
)
d = get_vertex_index(
i + 1,
j + 1,
k
)
face_append(
(a, b, c, d)
)
# ========================================================
# -Z
# ========================================================
if (
k == 0 or
not contains(
key - voxel_layer
)
):
a = get_vertex_index(
i,
j,
k
)
b = get_vertex_index(
i,
j + 1,
k
)
c = get_vertex_index(
i + 1,
j + 1,
k
)
d = get_vertex_index(
i + 1,
j,
k
)
face_append(
(a, b, c, d)
)
# ========================================================
# +Z
# ========================================================
if (
k == nz - 1 or
not contains(
key + voxel_layer
)
):
a = get_vertex_index(
i,
j,
k + 1
)
b = get_vertex_index(
i + 1,
j,
k + 1
)
c = get_vertex_index(
i + 1,
j + 1,
k + 1
)
d = get_vertex_index(
i,
j + 1,
k + 1
)
face_append(
(a, b, c, d)
)
boundary_time = (
time.time() -
boundary_start
)
print(
"Boundary scan : {:.2f} sec".format(
boundary_time
)
)
print(
"Vertices : {:,}".format(
len(vertex_keys)
)
)
print(
"Exposed faces : {:,}".format(
len(face_indices)
)
)
# vertex_map no longer needed
vertex_map = None
# ============================================================
# CREATE Point3f ARRAY
# ============================================================
print("")
print("Preparing compact vertex array...")
vertex_start = time.time()
points = List[
Rhino.Geometry.Point3f
]()
points.Capacity = len(
vertex_keys
)
point_add = points.Add
# Cache locals
sx = float(voxel_size)
ox = float(xmin)
oy = float(ymin)
oz = float(zmin)
# Decode grid keys directly
for key in vertex_keys:
iz = key // gv_layer
rem = key - iz * gv_layer
iy = rem // gvx
ix = rem - iy * gvx
point_add(
Rhino.Geometry.Point3f(
float(
ox +
ix * sx
),
float(
oy +
iy * sx
),
float(
oz +
iz * sx
)
)
)
vertex_keys = None
vertex_time = (
time.time() -
vertex_start
)
print(
"Vertex array : {:.2f} sec".format(
vertex_time
)
)
# ============================================================
# CREATE MeshFace ARRAY
# ============================================================
print("")
print("Preparing face array...")
face_start = time.time()
mesh_faces = List[
Rhino.Geometry.MeshFace
]()
mesh_faces.Capacity = len(
face_indices
)
face_add = mesh_faces.Add
for a, b, c, d in face_indices:
face_add(
Rhino.Geometry.MeshFace(
a,
b,
c,
d
)
)
face_indices = None
face_time = (
time.time() -
face_start
)
print(
"Face array : {:.2f} sec".format(
face_time
)
)
# ============================================================
# BULK CREATE RHINO MESH
# ============================================================
print("")
print("Creating Rhino mesh...")
rhino_start = time.time()
result = Rhino.Geometry.Mesh()
result.Vertices.Capacity = (
points.Count
)
result.Faces.Capacity = (
mesh_faces.Count
)
result.Vertices.AddVertices(
points
)
result.Faces.AddFaces(
mesh_faces
)
points = None
mesh_faces = None
result.FaceNormals.ComputeFaceNormals()
result.Normals.ComputeNormals()
result.Compact()
rhino_time = (
time.time() -
rhino_start
)
print(
"Rhino mesh : {:.2f} sec".format(
rhino_time
)
)
print(
"Output building : {:.2f} sec".format(
time.time() -
total_start
)
)
return result
# ================================================================
# MAIN
# ================================================================
def main():
ids = get_mesh_ids()
if not ids:
return
voxel_size = rs.GetReal(
"VoxelSize",
DEFAULT_VOXEL_SIZE,
0.000001
)
if voxel_size is None:
return
start_time = time.time()
old_redraw = rs.EnableRedraw(False)
try:
print("")
print("FAST VOXEL RASTERIZER")
print("=====================")
print(
"Meshes : {}".format(
len(ids)
)
)
print(
"VoxelSize : {}".format(
voxel_size
)
)
print(
"Workers : {}".format(
WORKERS
)
)
print(
"Sample step : {:.3f}".format(
voxel_size *
SAMPLE_FACTOR
)
)
# ========================================================
# COMBINE
# ========================================================
print("")
print("Combining meshes...")
mesh = combine_meshes(
ids
)
if mesh.Faces.Count == 0:
print(
"No valid mesh faces."
)
return
print(
"Triangles: {:,}".format(
mesh.Faces.Count
)
)
# ========================================================
# BOUNDING BOX
# ========================================================
bbox = mesh.GetBoundingBox(
True
)
padding = (
voxel_size *
PADDING_VOXELS
)
xmin = (
bbox.Min.X -
padding
)
ymin = (
bbox.Min.Y -
padding
)
zmin = (
bbox.Min.Z -
padding
)
xmax = (
bbox.Max.X +
padding
)
ymax = (
bbox.Max.Y +
padding
)
zmax = (
bbox.Max.Z +
padding
)
nx = max(
1,
int(
math.ceil(
(
xmax -
xmin
) /
voxel_size
)
)
)
ny = max(
1,
int(
math.ceil(
(
ymax -
ymin
) /
voxel_size
)
)
)
nz = max(
1,
int(
math.ceil(
(
zmax -
zmin
) /
voxel_size
)
)
)
print("")
print(
"Grid: {} x {} x {}".format(
nx,
ny,
nz
)
)
print(
"Full grid would contain {:,} voxels".format(
nx *
ny *
nz
)
)
# ========================================================
# EXTRACT TRIANGLES
# ========================================================
print("")
print(
"Extracting triangle coordinates..."
)
triangles = extract_triangles(
mesh
)
mesh = None
# ========================================================
# RASTERIZATION
# ========================================================
surface_voxels = rasterize_triangles(
triangles,
voxel_size,
xmin,
ymin,
zmin,
nx,
ny,
nz
)
triangles = None
print("")
print(
"Initial surface voxels: {:,}".format(
len(
surface_voxels
)
)
)
if not surface_voxels:
print(
"No voxels generated."
)
return
# ========================================================
# DILATION
# ========================================================
occupied = dilate_voxels(
surface_voxels,
nx,
ny,
nz,
DILATION_STEPS
)
surface_voxels = None
occupied_count = len(
occupied
)
print("")
print(
"Final occupied voxels: {:,}".format(
occupied_count
)
)
# ========================================================
# OUTPUT
# ========================================================
output = build_voxel_mesh_fast(
occupied,
xmin,
ymin,
zmin,
voxel_size,
nx,
ny,
nz
)
occupied = None
# ========================================================
# ADD TO RHINO
# ========================================================
print("")
print(
"Adding result to Rhino..."
)
result_id = sc.doc.Objects.AddMesh(
output
)
if (
result_id ==
System.Guid.Empty
):
print(
"Could not add output mesh."
)
return
elapsed = (
time.time() -
start_time
)
# ========================================================
# REPORT
# ========================================================
print("")
print("DONE")
print("====")
print(
"Occupied voxels : {:,}".format(
occupied_count
)
)
print(
"Vertices : {:,}".format(
output.Vertices.Count
)
)
print(
"Faces : {:,}".format(
output.Faces.Count
)
)
print(
"Closed : {}".format(
output.IsClosed
)
)
print(
"Valid : {}".format(
output.IsValid
)
)
print(
"TOTAL TIME : {:.2f} sec".format(
elapsed
)
)
finally:
rs.EnableRedraw(
old_redraw
)
sc.doc.Views.Redraw()
if __name__ == "__main__":
main()