Coverage for trimesh/voxel/creation.py: 71%
89 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1import numpy as np
3from .. import grouping, remesh, util
4from .. import transformations as tr
5from ..constants import log_time
6from ..typed import ArrayLike, Integer, Number, VoxelizationMethodsType
7from . import base
8from . import encoding as enc
11@log_time
12def voxelize_subdivide(
13 mesh, pitch: Number, max_iter: Integer | None = 10, edge_factor: Number = 2.0
14) -> base.VoxelGrid:
15 """
16 Voxelize a surface by subdividing a mesh until every edge is
17 shorter than: (pitch / edge_factor)
19 Parameters
20 -----------
21 mesh : trimesh.Trimesh
22 Source mesh
23 pitch
24 Side length of a single voxel cube
25 max_iter
26 Cap maximum subdivisions or None for no limit.
27 edge_factor
28 Proportion of pitch maximum edge length.
30 Returns
31 -----------
32 VoxelGrid instance representing the voxelized mesh.
33 """
34 max_edge = pitch / edge_factor
36 if max_iter is None:
37 longest_edge = np.linalg.norm(
38 mesh.vertices[mesh.edges[:, 0]] - mesh.vertices[mesh.edges[:, 1]], axis=1
39 ).max()
40 max_iter = max(int(np.ceil(np.log2(longest_edge / max_edge))), 0)
42 # get the same mesh sudivided so every edge is shorter
43 # than a factor of our pitch
44 v, _f, _idx = remesh.subdivide_to_size(
45 mesh.vertices, mesh.faces, max_edge=max_edge, max_iter=max_iter, return_index=True
46 )
48 # convert the vertices to their voxel grid position
49 # Provided edge_factor > 1 and max_iter is large enough, this is
50 # sufficient to preserve 6-connectivity at the level of voxels.
51 hit = np.round(v / pitch).astype(int)
53 # remove duplicates
54 unique, _inverse = grouping.unique_rows(hit)
56 # get the voxel centers in model space
57 occupied_index = hit[unique]
59 origin_index = occupied_index.min(axis=0)
60 origin_position = origin_index * pitch
62 return base.VoxelGrid(
63 enc.SparseBinaryEncoding(occupied_index - origin_index),
64 transform=tr.scale_and_translate(scale=pitch, translate=origin_position),
65 )
68def local_voxelize(
69 mesh,
70 point: ArrayLike,
71 pitch: Number,
72 radius: Number,
73 fill: bool = True,
74 **kwargs,
75) -> base.VoxelGrid | None:
76 """
77 Voxelize a mesh in the region of a cube around a point. When fill=True,
78 uses proximity.contains to fill the resulting voxels so may be meaningless
79 for non-watertight meshes. Useful to reduce memory cost for small values of
80 pitch as opposed to global voxelization.
82 Parameters
83 -----------
84 mesh : trimesh.Trimesh
85 Source geometry
86 point : (3, ) float
87 Point in space to voxelize around
88 pitch
89 Side length of a single voxel cube
90 radius
91 Number of voxel cubes to return in each direction.
92 kwargs
93 Parameters to pass to voxelize_subdivide
95 Returns
96 -----------
97 voxels : VoxelGrid instance with resolution (m, m, m) where m=2*radius+1
98 or None if the volume is empty
99 """
100 from scipy import ndimage
102 # make sure point is correct type/shape
103 point = np.asanyarray(point, dtype=np.float64).reshape(3)
104 # this is a gotcha- radius sounds a lot like it should be in
105 # float model space, not int voxel space so check
106 if not isinstance(radius, int):
107 raise ValueError("radius needs to be an integer number of cubes!")
109 # Bounds of region
110 bounds = np.concatenate(
111 (point - (radius + 0.5) * pitch, point + (radius + 0.5) * pitch)
112 )
114 # faces that intersect axis aligned bounding box
115 faces = list(mesh.triangles_tree.intersection(bounds))
117 # didn't hit anything so exit
118 if len(faces) == 0:
119 return None
121 local = mesh.submesh([[f] for f in faces], append=True)
123 # Translate mesh so point is at 0,0,0
124 local.apply_translation(-point)
126 # sparse, origin = voxelize_subdivide(local, pitch, **kwargs)
127 vox = voxelize_subdivide(local, pitch, **kwargs)
128 origin = vox.transform[:3, 3]
129 matrix = vox.encoding.dense
131 # Find voxel index for point
132 center = np.round(-origin / pitch).astype(np.int64)
134 # pad matrix if necessary
135 prepad = np.maximum(radius - center, 0)
136 postpad = np.maximum(center + radius + 1 - matrix.shape, 0)
138 matrix = np.pad(matrix, np.stack((prepad, postpad), axis=-1), mode="constant")
139 center += prepad
141 # Extract voxels within the bounding box
142 voxels = matrix[
143 center[0] - radius : center[0] + radius + 1,
144 center[1] - radius : center[1] + radius + 1,
145 center[2] - radius : center[2] + radius + 1,
146 ]
147 local_origin = point - radius * pitch # origin of local voxels
149 # Fill internal regions
150 if fill:
151 regions, n = ndimage.label(~voxels)
152 distance = ndimage.distance_transform_cdt(~voxels)
153 representatives = [
154 np.unravel_index((distance * (regions == i)).argmax(), distance.shape)
155 for i in range(1, n + 1)
156 ]
157 contains = mesh.contains(np.asarray(representatives) * pitch + local_origin)
158 where = np.where(contains)[0] + 1
159 internal = np.isin(regions.flatten(), where).reshape(regions.shape)
160 voxels = np.logical_or(voxels, internal)
162 return base.VoxelGrid(voxels, tr.translation_matrix(local_origin))
165@log_time
166def voxelize_ray(
167 mesh, pitch: Number, per_cell: ArrayLike | None = None
168) -> base.VoxelGrid:
169 """
170 Voxelize a mesh using ray queries.
172 Parameters
173 -------------
174 mesh
175 Mesh to be voxelized
176 pitch
177 Length of voxel cube
178 per_cell : (2,) int
179 How many ray queries to make per cell
181 Returns
182 -------------
183 grid
184 VoxelGrid instance representing the voxelized mesh.
185 """
186 if per_cell is None:
187 # how many rays per cell
188 per_cell = np.array([2, 2], dtype=np.int64)
189 else:
190 per_cell = np.array(per_cell, dtype=np.int64).reshape(2)
192 # edge length of cube voxels
193 pitch = float(pitch)
195 # create the ray origins in a grid
196 bounds = mesh.bounds[:, :2].copy()
197 # offset start so we get the requested number per cell
198 bounds[0] += pitch / (1.0 + per_cell)
199 # offset end so arange doesn't short us
200 bounds[1] += pitch
201 # on X we are doing multiple rays per voxel step
202 step = pitch / per_cell
203 # 2D grid
204 ray_ori = util.grid_arange(bounds, step=step)
205 # a Z position below the mesh
206 z = np.ones(len(ray_ori)) * (mesh.bounds[0][2] - pitch)
207 ray_ori = np.column_stack((ray_ori, z))
208 # all rays are along positive Z
209 ray_dir = np.ones_like(ray_ori) * [0, 0, 1]
211 # if you have pyembree this should be decently fast
212 hits = mesh.ray.intersects_location(ray_ori, ray_dir)[0]
214 # just convert hit locations to integer positions
215 voxels = np.round(hits / pitch).astype(np.int64)
217 # offset voxels by min, so matrix isn't huge
218 origin_index = voxels.min(axis=0)
219 voxels -= origin_index
220 encoding = enc.SparseBinaryEncoding(voxels)
221 origin_position = origin_index * pitch
222 return base.VoxelGrid(
223 encoding, tr.scale_and_translate(scale=pitch, translate=origin_position)
224 )
227@log_time
228def voxelize_binvox(
229 mesh,
230 pitch: Number | None = None,
231 dimension: Integer | None = None,
232 bounds: ArrayLike | None = None,
233 **binvoxer_kwargs,
234) -> base.VoxelGrid:
235 """
236 Voxelize via binvox tool.
238 Parameters
239 --------------
240 mesh : trimesh.Trimesh
241 Mesh to voxelize
242 pitch : float
243 Side length of each voxel. Ignored if dimension is provided
244 dimension: int
245 Number of voxels along each dimension. If not provided, this is
246 calculated based on pitch and bounds/mesh extents
247 bounds: (2, 3) float
248 min/max values of the returned `VoxelGrid` in each instance. Uses
249 `mesh.bounds` if not provided.
250 **binvoxer_kwargs:
251 Passed to `trimesh.exchange.binvox.Binvoxer`.
252 Should not contain `bounding_box` if bounds is not None.
254 Returns
255 --------------
256 grid
257 `VoxelGrid` instance
259 Raises
260 --------------
261 `ValueError` if `bounds is not None and 'bounding_box' in binvoxer_kwargs`.
262 """
263 from trimesh.exchange import binvox
265 if dimension is None:
266 # pitch must be provided
267 if bounds is None:
268 extents = mesh.extents
269 else:
270 mins, maxs = bounds
271 extents = maxs - mins
272 dimension = int(np.ceil(np.max(extents) / pitch))
273 if bounds is not None:
274 if "bounding_box" in binvoxer_kwargs:
275 raise ValueError("Cannot provide both bounds and bounding_box")
276 binvoxer_kwargs["bounding_box"] = np.asanyarray(bounds).flatten()
278 binvoxer = binvox.Binvoxer(dimension=dimension, **binvoxer_kwargs)
279 return binvox.voxelize_mesh(mesh, binvoxer)
282voxelizers = util.FunctionRegistry(
283 ray=voxelize_ray, subdivide=voxelize_subdivide, binvox=voxelize_binvox
284)
287def voxelize(
288 mesh,
289 pitch: Number | None,
290 method: VoxelizationMethodsType = "subdivide",
291 **kwargs,
292) -> base.VoxelGrid | None:
293 """
294 Voxelize the given mesh using the specified implementation.
296 See `voxelizers` for available implementations or to add your own, e.g. via
297 `voxelizers['custom_key'] = custom_fn`.
299 `custom_fn` should have signature `(mesh, pitch, **kwargs) -> VoxelGrid`
300 and should not modify encoding.
302 Parameters
303 --------------
304 mesh
305 Geometry to voxelize
306 pitch
307 Side length of each voxel.
308 method
309 Which voxelization method to use.
310 kwargs
311 Passed through to the specified implementation.
313 Returns
314 --------------
315 grid
316 A VoxelGrid instance.
317 """
318 return voxelizers(method, mesh=mesh, pitch=pitch, **kwargs)