Coverage for trimesh/bounds.py: 77%
199 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1import numpy as np
3from . import convex, geometry, grouping, nsphere, transformations, util
4from .constants import log, now
5from .typed import ArrayLike, NDArray
7try:
8 # scipy is a soft dependency
9 from scipy import optimize
10 from scipy.spatial import ConvexHull
11except BaseException as E:
12 # raise the exception when someone tries to use it
13 from . import exceptions
15 ConvexHull = exceptions.ExceptionWrapper(E)
16 optimize = exceptions.ExceptionWrapper(E)
18try:
19 from scipy.spatial import QhullError
20except BaseException:
21 QhullError = BaseException
23# a 90 degree rotation
24_flip = transformations.planar_matrix(theta=np.pi / 2)
25_flip.flags.writeable = False
28def oriented_bounds_2D(points, qhull_options="QbB"):
29 """
30 Find an oriented bounding box for an array of 2D points.
32 Details on qhull options:
33 http://www.qhull.org/html/qh-quick.htm#options
35 Parameters
36 ----------
37 points : (n,2) float
38 Points in 2D.
40 Returns
41 ----------
42 transform : (3,3) float
43 Homogeneous 2D transformation matrix to move the
44 input points so that the axis aligned bounding box
45 is CENTERED AT THE ORIGIN.
46 rectangle : (2,) float
47 Size of extents once input points are transformed
48 by transform
49 """
50 # create a convex hull object of our points
51 # 'QbB' is a qhull option which has it scale the input to unit
52 # box to avoid precision issues with very large/small meshes
53 convex = ConvexHull(points, qhull_options=qhull_options)
55 # (n,2,3) line segments
56 hull_edges = convex.points[convex.simplices]
57 # (n,2) points on the convex hull
58 hull_points = convex.points[convex.vertices]
60 # unit vector direction of the edges of the hull polygon
61 # filter out zero- magnitude edges via check_valid
62 edge_vectors = hull_edges[:, 1] - hull_edges[:, 0]
63 edge_norm = np.sqrt(np.dot(edge_vectors**2, [1, 1]))
64 edge_nonzero = edge_norm > 1e-10
65 edge_vectors = edge_vectors[edge_nonzero] / edge_norm[edge_nonzero].reshape((-1, 1))
67 # create a set of perpendicular vectors
68 perp_vectors = np.fliplr(edge_vectors) * [-1.0, 1.0]
70 # find the projection of every hull point on every edge vector
71 # this does create a potentially gigantic n^2 array in memory,
72 # and there is the 'rotating calipers' algorithm which avoids this
73 # however, we have reduced n with a convex hull and numpy dot products
74 # are extremely fast so in practice this usually ends up being fine
75 x = np.dot(edge_vectors, hull_points.T)
76 y = np.dot(perp_vectors, hull_points.T)
78 # reduce the projections to maximum and minimum per edge vector
79 bounds = np.column_stack((x.min(axis=1), y.min(axis=1), x.max(axis=1), y.max(axis=1)))
81 # calculate the extents and area for each edge vector pair
82 extents = np.diff(bounds.reshape((-1, 2, 2)), axis=1).reshape((-1, 2))
83 area = np.prod(extents, axis=1)
84 area_min = area.argmin()
86 # (2,) float of smallest rectangle size
87 rectangle = extents[area_min]
89 # find the (3,3) homogeneous transformation which moves the input
90 # points to have a bounding box centered at the origin
91 offset = -bounds[area_min][:2] - (rectangle * 0.5)
92 theta = np.arctan2(*edge_vectors[area_min][::-1])
93 transform = transformations.planar_matrix(offset, theta)
95 # we would like to consistently return an OBB with
96 # the largest dimension along the X axis rather than
97 # the long axis being arbitrarily X or Y.
98 if rectangle[1] > rectangle[0]:
99 # apply the rotation
100 transform = np.dot(_flip, transform)
101 # switch X and Y in the OBB extents
102 rectangle = rectangle[::-1]
104 return transform, rectangle
107def oriented_bounds(obj, angle_digits=1, ordered=True, normal=None, coplanar_tol=1e-12):
108 """
109 Find the oriented bounding box for a Trimesh
111 Parameters
112 ----------
113 obj : trimesh.Trimesh, (n, 2) float, or (n, 3) float
114 Mesh object or points in 2D or 3D space
115 angle_digits : int
116 How much angular precision do we want on our result.
117 Even with less precision the returned extents will cover
118 the mesh albeit with larger than minimal volume, and may
119 experience substantial speedups.
120 ordered : bool
121 Return a consistent order for bounds
122 normal : None or (3,) float
123 Override search for normal on 3D meshes.
124 coplanar_tol : float
125 If a convex hull fails and we are checking to see if the
126 points are coplanar this is the maximum deviation from
127 a plane where the points will be considered coplanar.
129 Returns
130 ----------
131 to_origin : (4,4) float
132 Transformation matrix which will move the center of the
133 bounding box of the input mesh to the origin.
134 extents: (3,) float
135 The extents of the mesh once transformed with to_origin
136 """
138 def oriented_bounds_coplanar(points):
139 """
140 Find an oriented bounding box for an array of coplanar 3D points.
142 Parameters
143 ----------
144 points : (n, 3) float
145 Points in 3D that occupy a 2D subspace.
147 Returns
148 ----------
149 to_origin : (4, 4) float
150 Transformation matrix which will move the center of the
151 bounding box of the input mesh to the origin.
152 extents : (3,) float
153 The extents of the mesh once transformed with to_origin
154 """
155 # Shift points about the origin and rotate into the xy plane
156 points_mean = np.mean(points, axis=0)
157 points_demeaned = points - points_mean
158 _, _, vh = np.linalg.svd(points_demeaned, full_matrices=False)
159 points_2d = np.matmul(points_demeaned, vh.T)
160 if np.any(np.abs(points_2d[:, 2]) > coplanar_tol):
161 raise ValueError("Points must be coplanar")
163 # Construct a homogeneous matrix representing the transformation above
164 to_2d = np.eye(4)
165 to_2d[:3, :3] = vh
166 to_2d[:3, 3] = -np.matmul(vh, points_mean)
168 # Find the 2D bounding box using the polygon
169 to_origin_2d, extents_2d = oriented_bounds_2D(points_2d[:, :2])
170 # Make extents 3D
171 extents = np.append(extents_2d, 0.0)
172 # convert transformation from 2D to 3D and combine
173 to_origin = np.matmul(transformations.planar_matrix_to_3D(to_origin_2d), to_2d)
174 return to_origin, extents
176 try:
177 # extract a set of convex hull vertices and normals from the input
178 # we bother to do this to avoid recomputing the full convex hull if
179 # possible
180 if hasattr(obj, "convex_hull"):
181 # if we have been passed a mesh, use its existing convex hull to pull from
182 # cache rather than recomputing. This version of the cached convex hull has
183 # normals pointing in arbitrary directions (straight from qhull)
184 # using this avoids having to compute the expensive corrected normals
185 # that mesh.convex_hull uses since normal directions don't matter
186 # here
187 hull = obj.convex_hull
188 elif util.is_sequence(obj):
189 # we've been passed a list of points
190 points = np.asanyarray(obj)
191 if util.is_shape(points, (-1, 2)):
192 return oriented_bounds_2D(points)
193 elif util.is_shape(points, (-1, 3)):
194 hull = convex.convex_hull(points, repair=True)
195 else:
196 raise ValueError("Points are not (n,3) or (n,2)!")
197 else:
198 raise ValueError("Oriented bounds must be passed a mesh or a set of points!")
199 except QhullError:
200 # Try to recover from Qhull error if due to mesh being less than 3
201 # dimensional
202 if hasattr(obj, "vertices"):
203 points = obj.vertices.view(np.ndarray)
204 elif util.is_sequence(obj):
205 points = np.asanyarray(obj)
206 else:
207 raise
208 return oriented_bounds_coplanar(points)
210 vertices = hull.vertices
211 hull_adj = hull.face_adjacency.T
212 hull_edge = hull.face_adjacency_edges
213 hull_normals = hull.face_normals
215 # matrices which will rotate each hull normal to [0,0,1]
216 if normal is None:
217 # convert face normals to spherical coordinates on the upper hemisphere
218 # the vector_hemisphere call effectively merges negative but otherwise
219 # identical vectors
220 spherical_coords = util.vector_to_spherical(util.vector_hemisphere(hull_normals))
221 # the unique_rows call on merge angles gets unique spherical directions to check
222 # we get a substantial speedup in the transformation matrix creation
223 # inside the loop by converting to angles ahead of time
224 spherical_unique = grouping.unique_rows(spherical_coords, digits=angle_digits)[0]
225 matrices = [
226 transformations.spherical_matrix(*s).T
227 for s in spherical_coords[spherical_unique]
228 ]
229 normals = util.spherical_to_vector(spherical_coords[spherical_unique])
230 else:
231 # if explicit normal was passed use it and skip the grouping
232 matrices = [geometry.align_vectors(normal, [0, 0, 1])]
233 normals = [normal]
235 tic = now()
236 min_2D = None
237 min_volume = np.inf
239 # we now need to loop through all the possible candidate
240 # directions for aligning our oriented bounding box.
241 for normal, to_2D in zip(normals, matrices):
242 # we could compute the hull in 2D for every direction
243 # but since we know we're dealing with a convex blob
244 # we can do back-face culling and then take the boundary
245 # start by picking the normal direction with fewer edges
246 side = np.dot(hull_normals, normal) > -1e-10
247 # for coplanar points this could be empty
248 if not side.any():
249 continue
250 # this line is a heavy lift as it is finding the pairs of
251 # adjacent faces where *exactly one* out of two of the faces
252 # is visible (xor) and then using the index to get the edge
253 edges = hull_edge[np.bitwise_xor(*side[hull_adj])]
255 # project the 3D convex hull vertices onto the plane
256 projected = np.dot(to_2D[:3, :3], vertices.T).T[:, :3]
257 # get the line segments of edges in 2D
258 edge_vert = projected[:, :2][edges]
259 # now get them as unit vectors
260 edge_vectors = edge_vert[:, 1, :] - edge_vert[:, 0, :]
261 edge_norm = np.sqrt(np.dot(edge_vectors**2, [1, 1]))
262 edge_nonzero = edge_norm > 1e-10
263 edge_vectors = edge_vectors[edge_nonzero] / edge_norm[edge_nonzero].reshape(
264 (-1, 1)
265 )
266 # create a set of perpendicular vectors
267 perp_vectors = np.fliplr(edge_vectors) * [-1.0, 1.0]
269 # find the projection of every hull point on every edge vector
270 # this does create a potentially gigantic n^2 array in memory
271 # and there is the 'rotating calipers' algorithm which avoids this
272 # however, we have reduced n with a convex hull and numpy dot products
273 # are extremely fast so in practice this usually ends up being fine
274 x = np.dot(edge_vectors, edge_vert[:, 0, :2].T)
275 y = np.dot(perp_vectors, edge_vert[:, 0, :2].T)
276 area = ((x.max(axis=1) - x.min(axis=1)) * (y.max(axis=1) - y.min(axis=1))).min()
278 # the volume is 2D area plus the projected height
279 volume = area * np.ptp(projected[:, 2])
281 # store this transform if it's better than one we've seen
282 if volume < min_volume:
283 min_volume = volume
284 min_2D = to_2D
286 # we know the minimum volume transform which should be the expensive
287 # part so now we need to do the bookkeeping to find the box
288 vert_ones = np.column_stack((vertices, np.ones(len(vertices)))).T
289 projected = np.dot(min_2D, vert_ones).T[:, :3]
290 height = np.ptp(projected[:, 2])
291 rotation_2D, box = oriented_bounds_2D(projected[:, :2])
292 min_extents = np.append(box, height)
293 rotation_2D[:2, 2] = 0.0
294 rotation_Z = transformations.planar_matrix_to_3D(rotation_2D)
296 # combine the 2D OBB transformation with the 2D projection transform
297 to_origin = np.dot(rotation_Z, min_2D)
299 # transform points using our matrix to find the translation
300 transformed = transformations.transform_points(vertices, to_origin)
301 box_center = transformed.min(axis=0) + np.ptp(transformed, axis=0) * 0.5
302 to_origin[:3, 3] = -box_center
304 # return ordered 3D extents
305 if ordered:
306 # sort the three extents
307 order = min_extents.argsort()
308 # generate a matrix which will flip transform
309 # to match the new ordering
310 flip = np.eye(4)
311 flip[:3, :3] = -np.eye(3)[order]
313 # make sure transform isn't mangling triangles
314 # by reversing windings on triangles
315 if not np.isclose(np.linalg.det(flip[:3, :3]), 1.0):
316 flip[:3, :3] = np.dot(flip[:3, :3], -np.eye(3))
318 # apply the flip to the OBB transform
319 to_origin = np.dot(flip, to_origin)
320 # apply the order to the extents
321 min_extents = min_extents[order]
323 log.debug("oriented_bounds checked %d vectors in %0.4fs", len(matrices), now() - tic)
325 return to_origin, min_extents
328def minimum_cylinder(obj, sample_count=6, angle_tol=0.001):
329 """
330 Find the approximate minimum volume cylinder which contains
331 a mesh or a a list of points.
333 Samples a hemisphere then uses scipy.optimize to pick the
334 final orientation of the cylinder.
336 A nice discussion about better ways to implement this is here:
337 https://www.staff.uni-mainz.de/schoemer/publications/ALGO00.pdf
340 Parameters
341 ----------
342 obj : trimesh.Trimesh, or (n, 3) float
343 Mesh object or points in space
344 sample_count : int
345 How densely should we sample the hemisphere.
346 Angular spacing is 180 degrees / this number
348 Returns
349 ----------
350 result : dict
351 With keys:
352 'radius' : float, radius of cylinder
353 'height' : float, height of cylinder
354 'transform' : (4,4) float, transform from the origin
355 to centered cylinder
356 """
358 def volume_from_angles(spherical, return_data=False):
359 """
360 Takes spherical coordinates and calculates the volume
361 of a cylinder along that vector
363 Parameters
364 ---------
365 spherical : (2,) float
366 Theta and phi
367 return_data : bool
368 Flag for returned
370 Returns
371 --------
372 if return_data:
373 transform ((4,4) float)
374 radius (float)
375 height (float)
376 else:
377 volume (float)
378 """
379 to_2D = transformations.spherical_matrix(*spherical, axes="rxyz")
380 projected = transformations.transform_points(hull, matrix=to_2D)
381 height = np.ptp(projected[:, 2])
383 try:
384 center_2D, radius = nsphere.minimum_nsphere(projected[:, :2])
385 except ValueError:
386 return np.inf
388 volume = np.pi * height * (radius**2)
389 if return_data:
390 center_3D = np.append(center_2D, projected[:, 2].min() + (height * 0.5))
391 transform = np.dot(
392 np.linalg.inv(to_2D), transformations.translation_matrix(center_3D)
393 )
394 return transform, radius, height
395 return volume
397 # we've been passed a mesh with radial symmetry
398 # use center mass and symmetry axis and go home early
399 if hasattr(obj, "symmetry") and obj.symmetry == "radial":
400 # find our origin
401 if obj.is_watertight:
402 # set origin to center of mass
403 origin = obj.center_mass
404 else:
405 # convex hull should be watertight
406 origin = obj.convex_hull.center_mass
407 # will align symmetry axis with Z and move origin to zero
408 to_2D = geometry.plane_transform(origin=origin, normal=obj.symmetry_axis)
409 # transform vertices to plane to check
410 on_plane = transformations.transform_points(obj.vertices, to_2D)
411 # cylinder height is overall Z span
412 height = np.ptp(on_plane[:, 2])
413 # center mass is correct on plane, but position
414 # along symmetry axis may be wrong so slide it
415 slide = transformations.translation_matrix(
416 [0, 0, (height / 2.0) - on_plane[:, 2].max()]
417 )
418 to_2D = np.dot(slide, to_2D)
419 # radius is maximum radius
420 radius = (on_plane[:, :2] ** 2).sum(axis=1).max() ** 0.5
421 # save kwargs
422 result = {"height": height, "radius": radius, "transform": np.linalg.inv(to_2D)}
423 return result
425 # get the points on the convex hull of the result
426 hull = convex.hull_points(obj)
427 if not util.is_shape(hull, (-1, 3)):
428 raise ValueError("Input must be reducable to 3D points!")
430 # sample a hemisphere so local hill climbing can do its thing
431 samples = util.grid_linspace([[0, 0], [np.pi, np.pi]], sample_count)
433 # if it's rotationally symmetric the bounding cylinder
434 # is almost certainly along one of the PCI vectors
435 if hasattr(obj, "principal_inertia_vectors"):
436 # add the principal inertia vectors if we have a mesh
437 samples = np.vstack(
438 (samples, util.vector_to_spherical(obj.principal_inertia_vectors))
439 )
441 tic = [now()]
442 # the projected volume at each sample
443 volumes = np.array([volume_from_angles(i) for i in samples])
444 # the best vector in (2,) spherical coordinates
445 best = samples[volumes.argmin()]
446 tic.append(now())
448 # since we already explored the global space, set the bounds to be
449 # just around the sample that had the lowest volume
450 step = 2 * np.pi / sample_count
451 bounds = [(best[0] - step, best[0] + step), (best[1] - step, best[1] + step)]
452 # run the local optimization
453 r = optimize.minimize(
454 volume_from_angles, best, tol=angle_tol, method="SLSQP", bounds=bounds
455 )
457 tic.append(now())
458 log.debug("Performed search in %f and minimize in %f", *np.diff(tic))
460 # actually chunk the information about the cylinder
461 transform, radius, height = volume_from_angles(r["x"], return_data=True)
463 result = {"transform": transform, "radius": radius, "height": height}
464 return result
467def to_extents(bounds):
468 """
469 Convert an axis aligned bounding box to extents and
470 transform.
472 Parameters
473 ------------
474 bounds : (2, 3) float
475 Axis aligned bounds in space
477 Returns
478 ------------
479 extents : (3,) float
480 Extents of the bounding box
481 transform : (4, 4) float
482 Homogeneous transform moving extents to bounds
483 """
484 bounds = np.asanyarray(bounds, dtype=np.float64)
485 if bounds.shape != (2, 3):
486 raise ValueError("bounds must be (2, 3)")
488 extents = np.ptp(bounds, axis=0)
489 transform = np.eye(4)
490 transform[:3, 3] = bounds.mean(axis=0)
492 return extents, transform
495def corners(bounds):
496 """
497 Given a pair of axis aligned bounds, return all
498 8 corners of the bounding box.
500 Parameters
501 ----------
502 bounds : (2,3) or (2,2) float
503 Axis aligned bounds
505 Returns
506 ----------
507 corners : (8,3) float
508 Corner vertices of the cube
509 """
511 bounds = np.asanyarray(bounds, dtype=np.float64)
513 if util.is_shape(bounds, (2, 2)):
514 bounds = np.column_stack((bounds, [0, 0]))
515 elif not util.is_shape(bounds, (2, 3)):
516 raise ValueError("bounds must be (2,2) or (2,3)!")
518 minx, miny, minz, maxx, maxy, maxz = np.arange(6)
519 corner_index = np.array(
520 [
521 minx,
522 miny,
523 minz,
524 maxx,
525 miny,
526 minz,
527 maxx,
528 maxy,
529 minz,
530 minx,
531 maxy,
532 minz,
533 minx,
534 miny,
535 maxz,
536 maxx,
537 miny,
538 maxz,
539 maxx,
540 maxy,
541 maxz,
542 minx,
543 maxy,
544 maxz,
545 ]
546 ).reshape((-1, 3))
548 corners = bounds.reshape(-1)[corner_index]
549 return corners
552def contains(bounds: ArrayLike, points: ArrayLike) -> NDArray[np.bool_]:
553 """
554 Do an axis aligned bounding box check on an array of points.
556 Parameters
557 -----------
558 bounds : (2, dimension) float
559 Axis aligned bounding box
560 points : (n, dimension) float
561 Points in space
563 Returns
564 -----------
565 points_inside : (n,) bool
566 True if points are inside the AABB
567 """
568 # make sure we have correct input types
569 bounds = np.asanyarray(bounds, dtype=np.float64)
570 points = np.asanyarray(points, dtype=np.float64)
572 if len(bounds) != 2:
573 raise ValueError("bounds must be (2,dimension)!")
574 if not util.is_shape(points, (-1, bounds.shape[1])):
575 raise ValueError("bounds shape must match points!")
577 # run the simple check
578 points_inside = np.logical_and(
579 (points > bounds[0]).all(axis=1), (points < bounds[1]).all(axis=1)
580 )
582 return points_inside