Coverage for trimesh/voxel/ops.py: 78%
186 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1import numpy as np
3from .. import util
4from ..constants import log
5from ..typed import ArrayLike, Number
8def fill_orthographic(dense):
9 shape = dense.shape
10 indices = np.stack(
11 np.meshgrid(*(np.arange(s) for s in shape), indexing="ij"), axis=-1
12 )
13 empty = np.logical_not(dense)
15 def fill_axis(axis):
16 base_local_indices = indices[..., axis]
17 local_indices = base_local_indices.copy()
18 local_indices[empty] = shape[axis]
19 mins = np.min(local_indices, axis=axis, keepdims=True)
20 local_indices = base_local_indices.copy()
21 local_indices[empty] = -1
22 maxs = np.max(local_indices, axis=axis, keepdims=True)
24 return np.logical_and(
25 base_local_indices >= mins,
26 base_local_indices <= maxs,
27 )
29 filled = fill_axis(axis=0)
30 for axis in range(1, len(shape)):
31 filled = np.logical_and(filled, fill_axis(axis))
32 return filled
35def fill_base(sparse_indices):
36 """
37 Given a sparse surface voxelization, fill in between columns.
39 Parameters
40 --------------
41 sparse_indices: (n, 3) int, location of filled cells
43 Returns
44 --------------
45 filled: (m, 3) int, location of filled cells
46 """
47 # validate inputs
48 sparse_indices = np.asanyarray(sparse_indices, dtype=np.int64)
49 if not util.is_shape(sparse_indices, (-1, 3)):
50 raise ValueError("incorrect shape")
52 # create grid and mark inner voxels
53 max_value = sparse_indices.max() + 3
55 grid = np.zeros((max_value, max_value, max_value), bool)
56 voxels_sparse = np.add(sparse_indices, 1)
58 grid[tuple(voxels_sparse.T)] = 1
60 for i in range(max_value):
61 check_dir2 = False
62 for j in range(0, max_value - 1):
63 idx = []
64 # find transitions first
65 # transition positions are from 0 to 1 and from 1 to 0
66 eq = np.equal(grid[i, j, :-1], grid[i, j, 1:])
67 idx = np.where(np.logical_not(eq))[0] + 1
68 c = len(idx)
69 check_dir2 = (c % 4) > 0 and c > 4
70 if c < 4:
71 continue
72 for s in range(0, c - c % 4, 4):
73 grid[i, j, idx[s] : idx[s + 3]] = 1
74 if not check_dir2:
75 continue
77 # check another direction for robustness
78 for k in range(0, max_value - 1):
79 idx = []
80 # find transitions first
81 eq = np.equal(grid[i, :-1, k], grid[i, 1:, k])
82 idx = np.where(np.logical_not(eq))[0] + 1
83 c = len(idx)
84 if c < 4:
85 continue
86 for s in range(0, c - c % 4, 4):
87 grid[i, idx[s] : idx[s + 3], k] = 1
89 # generate new voxels
90 filled = np.column_stack(np.where(grid))
91 filled -= 1
93 return filled
96fill_voxelization = fill_base
99def matrix_to_marching_cubes(
100 matrix: ArrayLike,
101 pitch: Number | ArrayLike = 1.0,
102 threshold: Number | None = None,
103):
104 """
105 Convert an (n, m, p) matrix into a mesh, using marching_cubes.
107 Parameters
108 -----------
109 matrix : (n, m, p) bool
110 Occupancy array
111 pitch : float or length-3 tuple of floats, optional
112 Voxel spacing in each dimension
113 threshold : float or None, optional
114 If specified, converts the input into a boolean
115 matrix by considering values above `threshold` as True
118 Returns
119 ----------
120 mesh : trimesh.Trimesh
121 Mesh generated by meshing voxels using
122 the marching cubes algorithm in skimage
123 """
124 from skimage import measure
126 from ..base import Trimesh
128 if threshold is not None:
129 matrix = np.asanyarray(matrix) > threshold
130 else:
131 matrix = np.asanyarray(matrix, dtype=bool)
133 rev_matrix = np.logical_not(matrix) # Takes set about 0.
134 # Add in padding so marching cubes can function properly with
135 # voxels on edge of AABB
136 pad_width = 1
137 rev_matrix = np.pad(
138 rev_matrix, pad_width=(pad_width), mode="constant", constant_values=(1)
139 )
141 # pick between old and new API
142 if hasattr(measure, "marching_cubes_lewiner"):
143 func = measure.marching_cubes_lewiner
144 else:
145 func = measure.marching_cubes
147 # Run marching cubes.
148 pitch = np.asanyarray(pitch)
149 if pitch.size == 1:
150 pitch = (pitch,) * 3
151 meshed = func(
152 volume=rev_matrix,
153 level=0.5,
154 spacing=pitch, # it is a boolean voxel grid
155 )
157 # allow results from either marching cubes function in skimage
158 # binaries available for python 3.3 and 3.4 appear to use the classic
159 # method
160 if len(meshed) == 2:
161 log.warning("using old marching cubes, may not be watertight!")
162 vertices, faces = meshed
163 normals = None
164 elif len(meshed) == 4:
165 vertices, faces, normals, _vals = meshed
167 # Return to the origin, add in the pad_width
168 vertices = np.subtract(vertices, pad_width * pitch)
169 # create the mesh
170 mesh = Trimesh(vertices=vertices, faces=faces, vertex_normals=normals)
171 return mesh
174def sparse_to_matrix(sparse):
175 """
176 Take a sparse (n,3) list of integer indexes of filled cells,
177 turn it into a dense (m,o,p) matrix.
179 Parameters
180 -----------
181 sparse : (n, 3) int
182 Index of filled cells
184 Returns
185 ------------
186 dense : (m, o, p) bool
187 Matrix of filled cells
188 """
190 sparse = np.asanyarray(sparse, dtype=np.int64)
191 if not util.is_shape(sparse, (-1, 3)):
192 raise ValueError("sparse must be (n,3)!")
194 shape = sparse.max(axis=0) + 1
195 matrix = np.zeros(np.prod(shape), dtype=bool)
196 multiplier = np.array([np.prod(shape[1:]), shape[2], 1])
198 index = (sparse * multiplier).sum(axis=1)
199 matrix[index] = True
201 dense = matrix.reshape(shape)
202 return dense
205def points_to_marching_cubes(points, pitch=1.0):
206 """
207 Mesh points by assuming they fill a voxel box, and then
208 running marching cubes on them
210 Parameters
211 ------------
212 points : (n, 3) float
213 Points in 3D space
215 Returns
216 -------------
217 mesh : trimesh.Trimesh
218 Points meshed using marching cubes
219 """
220 # make sure inputs are as expected
221 points = np.asanyarray(points, dtype=np.float64)
222 pitch = np.asanyarray(pitch, dtype=float)
224 # find the minimum value of points for origin
225 origin = points.min(axis=0)
226 # convert points to occupied voxel cells
227 index = ((points - origin) / pitch).round().astype(np.int64)
229 # convert voxel indices to a matrix
230 matrix = sparse_to_matrix(index)
232 # run marching cubes on the matrix to generate a mesh
233 mesh = matrix_to_marching_cubes(matrix, pitch=pitch)
234 mesh.vertices += origin
236 return mesh
239def multibox(centers, pitch=1.0, colors=None, remove_internal_faces=False):
240 """
241 Return a Trimesh object with a box at every center.
243 Doesn't do anything nice or fancy.
245 Parameters
246 -----------
247 centers : (n, 3) float
248 Center of boxes that are occupied
249 pitch : float
250 The edge length of a voxel
251 colors : (3,) or (4,) or (n,3) or (n, 4) float
252 Color of boxes
253 remove_internal_faces : bool
254 If True, removes internal faces shared between adjacent boxes
256 Returns
257 ---------
258 rough : Trimesh
259 Mesh object representing inputs
260 """
261 from .. import primitives
262 from ..base import Trimesh
264 # get centers as numpy array
265 centers = np.asanyarray(centers, dtype=np.float64)
267 # get a basic box
268 b = primitives.Box()
269 # apply the pitch
270 b.apply_scale(float(pitch))
271 # tile into one box vertex per center
272 v = np.tile(centers, (1, len(b.vertices))).reshape((-1, 3))
273 # offset to centers
274 v += np.tile(b.vertices, (len(centers), 1))
276 f = np.tile(b.faces, (len(centers), 1))
277 f += np.repeat(np.arange(len(centers)) * len(b.vertices), len(b.faces))[:, None]
279 if remove_internal_faces:
280 # Get 12 unit normals (1 per triangle face) indicating face direction
281 base_normals = np.round(b.face_normals).astype(int) # (12, 3)
282 # Expand those directions across all voxel boxes so as to check neighbor presence
283 face_normals = np.tile(base_normals, (len(centers), 1)) # (len(centers) * 12, 3)
285 # Maps each face to the voxel box it came from
286 face_voxel_idx = np.repeat(
287 np.arange(len(centers)), len(b.faces)
288 ) # (len(centers) * 12, )
289 # Converts voxel centers to discrete grid coordinates
290 voxel_coords = np.round(centers / pitch).astype(int) # (len(centers), 3)
291 # Creates a fast lookup structure for checking voxel neighbors
292 voxel_set = set(map(tuple, voxel_coords))
294 # Gets the grid coordinate of the voxel that owns each face
295 voxel_face_coords = voxel_coords[face_voxel_idx]
296 # Computes the adjacent voxel coordinate in the face direction
297 neighbor_coords = voxel_face_coords + face_normals
298 # Keeps only faces whose neighboring voxel does not exist
299 keep_mask = np.array([tuple(c) not in voxel_set for c in neighbor_coords])
300 else:
301 keep_mask = np.ones(len(f), dtype=bool)
303 face_colors = None
304 if colors is not None:
305 colors = np.asarray(colors)
306 if colors.ndim == 1:
307 colors = colors[None].repeat(len(centers), axis=0)
308 if colors.ndim == 2 and len(colors) == len(centers):
309 face_colors = colors.repeat(12, axis=0)[keep_mask]
311 mesh = Trimesh(vertices=v, faces=f[keep_mask], face_colors=face_colors)
313 return mesh
316def boolean_sparse(a, b, operation=np.logical_and):
317 """
318 Find common rows between two arrays very quickly
319 using 3D boolean sparse matrices.
321 Parameters
322 -----------
323 a: (n, d) int, coordinates in space
324 b: (m, d) int, coordinates in space
325 operation: numpy operation function, ie:
326 np.logical_and
327 np.logical_or
329 Returns
330 -----------
331 coords: (q, d) int, coordinates in space
332 """
333 # 3D sparse arrays, using wrapped scipy.sparse
334 # pip install sparse
335 import sparse
337 # find the bounding box of both arrays
338 extrema = np.array([a.min(axis=0), a.max(axis=0), b.min(axis=0), b.max(axis=0)])
339 origin = extrema.min(axis=0) - 1
340 size = tuple(np.ptp(extrema, axis=0) + 2)
342 # put nearby voxel arrays into same shape sparse array
343 sp_a = sparse.COO((a - origin).T, data=np.ones(len(a), dtype=bool), shape=size)
344 sp_b = sparse.COO((b - origin).T, data=np.ones(len(b), dtype=bool), shape=size)
346 # apply the logical operation
347 # get a sparse matrix out
348 applied = operation(sp_a, sp_b)
349 # reconstruct the original coordinates
350 coords = np.column_stack(applied.coords) + origin
352 return coords
355def strip_array(data):
356 shape = data.shape
357 ndims = len(shape)
358 padding = []
359 slices = []
360 for dim in range(len(shape)):
361 axis = tuple(range(dim)) + tuple(range(dim + 1, ndims))
362 filled = np.any(data, axis=axis)
363 (indices,) = np.nonzero(filled)
364 pad_left = indices[0]
365 pad_right = indices[-1]
366 padding.append([pad_left, pad_right])
367 slices.append(slice(pad_left, pad_right))
368 return data[tuple(slices)], np.array(padding, int)
371def indices_to_points(indices, pitch=None, origin=None):
372 """
373 Convert indices of an (n,m,p) matrix into a set of voxel center points.
375 Parameters
376 ----------
377 indices: (q, 3) int, index of voxel matrix (n,m,p)
378 pitch: float, what pitch was the voxel matrix computed with
379 origin: (3,) float, what is the origin of the voxel matrix
381 Returns
382 ----------
383 points: (q, 3) float, list of points
384 """
385 indices = np.asanyarray(indices)
386 if indices.shape[1:] != (3,):
387 raise ValueError("shape of indices must be (q, 3)")
389 points = np.array(indices, dtype=np.float64)
390 if pitch is not None:
391 points *= float(pitch)
392 if origin is not None:
393 origin = np.asanyarray(origin)
394 if origin.shape != (3,):
395 raise ValueError("shape of origin must be (3,)")
396 points += origin
398 return points
401def matrix_to_points(matrix, pitch=None, origin=None):
402 """
403 Convert an (n,m,p) matrix into a set of points for each voxel center.
405 Parameters
406 -----------
407 matrix: (n,m,p) bool, voxel matrix
408 pitch: float, what pitch was the voxel matrix computed with
409 origin: (3,) float, what is the origin of the voxel matrix
411 Returns
412 ----------
413 points: (q, 3) list of points
414 """
415 indices = np.column_stack(np.nonzero(matrix))
416 points = indices_to_points(indices=indices, pitch=pitch, origin=origin)
417 return points
420def points_to_indices(points, pitch=None, origin=None):
421 """
422 Convert center points of an (n,m,p) matrix into its indices.
424 Parameters
425 ----------
426 points : (q, 3) float
427 Center points of voxel matrix (n,m,p)
428 pitch : float
429 What pitch was the voxel matrix computed with
430 origin : (3,) float
431 What is the origin of the voxel matrix
433 Returns
434 ----------
435 indices : (q, 3) int
436 List of indices
437 """
438 points = np.array(points, dtype=np.float64)
439 if points.shape != (points.shape[0], 3):
440 raise ValueError("shape of points must be (q, 3)")
442 if origin is not None:
443 origin = np.asanyarray(origin)
444 if origin.shape != (3,):
445 raise ValueError("shape of origin must be (3,)")
446 points -= origin
447 if pitch is not None:
448 points /= pitch
450 origin = np.asanyarray(origin, dtype=np.float64)
451 pitch = float(pitch)
453 indices = np.round(points).astype(int)
454 return indices