Coverage for trimesh/intersections.py: 93%
258 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2intersections.py
3------------------
5Primarily mesh-plane intersections (slicing).
6"""
8import numpy as np
10from . import geometry, grouping, util
11from . import transformations as tf
12from . import triangles as tm
13from .constants import tol
14from .triangles import points_to_barycentric
17def mesh_plane(
18 mesh,
19 plane_normal,
20 plane_origin,
21 return_faces=False,
22 local_faces=None,
23 cached_dots=None,
24):
25 """
26 Find a the intersections between a mesh and a plane,
27 returning a set of line segments on that plane.
29 Parameters
30 ---------
31 mesh : Trimesh object
32 Source mesh to slice
33 plane_normal : (3,) float
34 Normal vector of plane to intersect with mesh
35 plane_origin : (3,) float
36 Point on plane to intersect with mesh
37 return_faces : bool
38 If True return face index each line is from
39 local_faces : None or (m,) int
40 Limit section to just these faces.
41 cached_dots : (n, 3) float
42 If an external function has stored dot
43 products pass them here to avoid recomputing.
45 Returns
46 ----------
47 lines : (m, 2, 3) float
48 List of 3D line segments in space.
49 face_index : (m,) int
50 Index of mesh.faces for each line
51 Only returned if return_faces was True
52 """
54 def triangle_cases(signs):
55 """
56 Figure out which faces correspond to which intersection
57 case from the signs of the dot product of each vertex.
58 Does this by bitbang each row of signs into an 8 bit
59 integer.
61 code : signs : intersects
62 0 : [-1 -1 -1] : No
63 2 : [-1 -1 0] : No
64 4 : [-1 -1 1] : Yes; 2 on one side, 1 on the other
65 6 : [-1 0 0] : Yes; one edge fully on plane
66 8 : [-1 0 1] : Yes; one vertex on plane 2 on different sides
67 12 : [-1 1 1] : Yes; 2 on one side, 1 on the other
68 14 : [0 0 0] : No (on plane fully)
69 16 : [0 0 1] : Yes; one edge fully on plane
70 20 : [0 1 1] : No
71 28 : [1 1 1] : No
73 Parameters
74 ----------
75 signs: (n,3) int, all values are -1,0, or 1
76 Each row contains the dot product of all three vertices
77 in a face with respect to the plane
79 Returns
80 ---------
81 basic : (n,) bool
82 Which faces are in the basic intersection case
83 one_vertex : (n,) bool
84 Which faces are in the one vertex case
85 one_edge : (n,) bool
86 Which faces are in the one edge case
87 """
89 signs_sorted = np.sort(signs, axis=1)
90 coded = np.zeros(len(signs_sorted), dtype=np.int8) + 14
91 for i in range(3):
92 coded += signs_sorted[:, i] << 3 - i
94 # one edge fully on the plane
95 # note that we are only accepting *one* of the on- edge cases,
96 # where the other vertex has a positive dot product (16) instead
97 # of both on- edge cases ([6, 16])
98 # this is so that for regions that are co-planar with the the section plane
99 # we don't end up with an invalid boundary
100 key = np.zeros(29, dtype=bool)
101 key[16] = True
102 one_edge = key[coded]
104 # one vertex on plane, other two on different sides
105 key[:] = False
106 key[8] = True
107 one_vertex = key[coded]
109 # one vertex on one side of the plane, two on the other
110 key[:] = False
111 key[[4, 12]] = True
112 basic = key[coded]
114 return basic, one_vertex, one_edge
116 def handle_on_vertex(signs, faces, vertices):
117 # case where one vertex is on plane
118 # and two are on different sides
119 vertex_plane = faces[signs == 0]
120 edge_thru = faces[signs != 0].reshape((-1, 2))
121 point_intersect, valid = plane_lines(
122 plane_origin, plane_normal, vertices[edge_thru.T], line_segments=False
123 )
124 lines = np.column_stack((vertices[vertex_plane[valid]], point_intersect)).reshape(
125 (-1, 2, 3)
126 )
127 return lines
129 def handle_on_edge(signs, faces, vertices):
130 # case where two vertices are on the plane and one is off
131 edges = faces[signs == 0].reshape((-1, 2))
132 points = vertices[edges]
133 return points
135 def handle_basic(signs, faces, vertices):
136 # case where one vertex is on one side and two are on the other
137 unique_element = grouping.unique_value_in_row(signs, unique=[-1, 1])
138 edges = np.column_stack(
139 (
140 faces[unique_element],
141 faces[np.roll(unique_element, 1, axis=1)],
142 faces[unique_element],
143 faces[np.roll(unique_element, 2, axis=1)],
144 )
145 ).reshape((-1, 2))
146 intersections, valid = plane_lines(
147 plane_origin, plane_normal, vertices[edges.T], line_segments=False
148 )
149 # since the data has been pre- culled, any invalid intersections at all
150 # means the culling was done incorrectly and thus things are broken
151 assert valid.all()
152 return intersections.reshape((-1, 2, 3))
154 # check input plane
155 plane_normal = np.asanyarray(plane_normal, dtype=np.float64)
156 plane_origin = np.asanyarray(plane_origin, dtype=np.float64)
157 if plane_origin.shape != (3,) or plane_normal.shape != (3,):
158 raise ValueError("Plane origin and normal must be (3,)!")
160 if local_faces is None:
161 # do a cross section against all faces
162 faces = mesh.faces
163 else:
164 local_faces = np.asanyarray(local_faces, dtype=np.int64)
165 # only take the subset of faces if passed
166 faces = mesh.faces[local_faces]
168 if cached_dots is not None:
169 dots = cached_dots
170 else:
171 # dot product of each vertex with the plane normal indexed by face
172 # so for each face the dot product of each vertex is a row
173 # shape is the same as mesh.faces (n,3)
174 dots = np.dot(mesh.vertices - plane_origin, plane_normal)
176 # sign of the dot product is -1, 0, or 1
177 # shape is the same as mesh.faces (n,3)
178 signs = np.zeros(len(mesh.vertices), dtype=np.int8)
179 signs[dots < -tol.merge] = -1
180 signs[dots > tol.merge] = 1
181 signs = signs[faces]
183 # figure out which triangles are in the cross section,
184 # and which of the three intersection cases they are in
185 cases = triangle_cases(signs)
186 # handlers for each case
187 handlers = (handle_basic, handle_on_vertex, handle_on_edge)
189 # the (m, 2, 3) line segments
190 lines = np.vstack(
191 [h(signs[c], faces[c], mesh.vertices) for c, h in zip(cases, handlers)]
192 )
194 if return_faces:
195 # everything that hit something
196 index = np.hstack([np.nonzero(c)[0] for c in cases])
197 assert index.dtype.kind == "i"
198 if local_faces is None:
199 return lines, index
200 # we are considering a subset of faces
201 # so we need to take the indexes from original
202 return lines, local_faces[index]
203 return lines
206def mesh_multiplane(mesh, plane_origin, plane_normal, heights):
207 """
208 A utility function for slicing a mesh by multiple
209 parallel planes which caches the dot product operation.
211 Parameters
212 -------------
213 mesh : trimesh.Trimesh
214 Geometry to be sliced by planes
215 plane_origin : (3,) float
216 Point on a plane
217 plane_normal : (3,) float
218 Normal vector of plane
219 heights : (m,) float
220 Offset distances from plane to slice at:
221 at `height=0` it will be exactly on the passed plane.
223 Returns
224 --------------
225 lines : (m,) sequence of (n, 2, 2) float
226 Lines in space for m planes
227 to_3D : (m, 4, 4) float
228 Transform to move each section back to 3D
229 face_index : (m,) sequence of (n,) int
230 Indexes of mesh.faces for each segment
231 """
232 # check input plane
233 plane_normal = util.unitize(plane_normal)
234 plane_origin = np.asanyarray(plane_origin, dtype=np.float64)
235 heights = np.asanyarray(heights, dtype=np.float64)
237 # dot product of every vertex with plane
238 vertex_dots = np.dot(plane_normal, (mesh.vertices - plane_origin).T)
240 # reconstruct transforms for each 2D section
241 base_transform = geometry.plane_transform(origin=plane_origin, normal=plane_normal)
242 base_transform = np.linalg.inv(base_transform)
244 # alter translation Z inside loop
245 translation = np.eye(4)
247 # store results
248 transforms = []
249 face_index = []
250 segments = []
252 # loop through user specified heights
253 for height in heights:
254 # offset the origin by the height
255 new_origin = plane_origin + (plane_normal * height)
256 # offset the dot products by height and index by faces
257 new_dots = vertex_dots - height
258 # run the intersection with the cached dot products
259 lines, index = mesh_plane(
260 mesh=mesh,
261 plane_origin=new_origin,
262 plane_normal=plane_normal,
263 return_faces=True,
264 cached_dots=new_dots,
265 )
267 # get the transforms to 3D space and back
268 translation[2, 3] = height
269 to_3D = np.dot(base_transform, translation)
270 to_2D = np.linalg.inv(to_3D)
271 transforms.append(to_3D)
273 # transform points to 2D frame
274 lines_2D = tf.transform_points(lines.reshape((-1, 3)), to_2D)
276 # if we didn't screw up the transform all
277 # of the Z values should be zero
278 # assert np.allclose(lines_2D[:, 2], 0.0)
280 # reshape back in to lines and discard Z
281 lines_2D = lines_2D[:, :2].reshape((-1, 2, 2))
282 # store (n, 2, 2) float lines
283 segments.append(lines_2D)
284 # store (n,) int indexes of mesh.faces
285 face_index.append(index)
287 # (n, 4, 4) transforms from 2D to 3D
288 transforms = np.array(transforms, dtype=np.float64)
290 return segments, transforms, face_index
293def plane_lines(plane_origin, plane_normal, endpoints, line_segments=True):
294 """
295 Calculate plane-line intersections
297 Parameters
298 ---------
299 plane_origin : (3,) float
300 Point on plane
301 plane_normal : (3,) float
302 Plane normal vector
303 endpoints : (2, n, 3) float
304 Points defining lines to be tested
305 line_segments : bool
306 If True, only returns intersections as valid if
307 vertices from endpoints are on different sides
308 of the plane.
310 Returns
311 ---------
312 intersections : (m, 3) float
313 Cartesian intersection points
314 valid : (n, 3) bool
315 Indicate whether a valid intersection exists
316 for each input line segment
317 """
318 endpoints = np.asanyarray(endpoints)
319 plane_origin = np.asanyarray(plane_origin).reshape(3)
320 line_dir = util.unitize(endpoints[1] - endpoints[0])
321 plane_normal = util.unitize(np.asanyarray(plane_normal).reshape(3))
323 t = np.dot(plane_normal, (plane_origin - endpoints[0]).T)
324 b = np.dot(plane_normal, line_dir.T)
326 # If the plane normal and line direction are perpendicular, it means
327 # the vector is 'on plane', and there isn't a valid intersection.
328 # We discard on-plane vectors by checking that the dot product is nonzero
329 valid = np.abs(b) > tol.zero
330 if line_segments:
331 test = np.dot(plane_normal, np.transpose(plane_origin - endpoints[1]))
332 different_sides = np.sign(t) != np.sign(test)
333 nonzero = np.logical_or(np.abs(t) > tol.zero, np.abs(test) > tol.zero)
334 valid = np.logical_and(valid, different_sides)
335 valid = np.logical_and(valid, nonzero)
337 d = np.divide(t[valid], b[valid])
338 intersection = endpoints[0][valid]
339 intersection = intersection + np.reshape(d, (-1, 1)) * line_dir[valid]
341 return intersection, valid
344def planes_lines(
345 plane_origins,
346 plane_normals,
347 line_origins,
348 line_directions,
349 return_distance=False,
350 return_denom=False,
351):
352 """
353 Given one line per plane find the intersection points.
355 Parameters
356 -----------
357 plane_origins : (n,3) float
358 Point on each plane
359 plane_normals : (n,3) float
360 Normal vector of each plane
361 line_origins : (n,3) float
362 Point at origin of each line
363 line_directions : (n,3) float
364 Direction vector of each line
365 return_distance : bool
366 Return distance from origin to point also
367 return_denom : bool
368 Return denominator, so you can check for small values
370 Returns
371 ----------
372 on_plane : (n,3) float
373 Points on specified planes
374 valid : (n,) bool
375 Did plane intersect line or not
376 distance : (n,) float
377 [OPTIONAL] Distance from point
378 denom : (n,) float
379 [OPTIONAL] Denominator
380 """
382 # check input types
383 plane_origins = np.asanyarray(plane_origins, dtype=np.float64)
384 plane_normals = np.asanyarray(plane_normals, dtype=np.float64)
385 line_origins = np.asanyarray(line_origins, dtype=np.float64)
386 line_directions = np.asanyarray(line_directions, dtype=np.float64)
388 # vector from line to plane
389 origin_vectors = plane_origins - line_origins
391 projection_ori = util.diagonal_dot(origin_vectors, plane_normals)
392 projection_dir = util.diagonal_dot(line_directions, plane_normals)
394 valid = np.abs(projection_dir) > 1e-5
396 distance = np.divide(projection_ori[valid], projection_dir[valid])
398 on_plane = line_directions[valid] * distance.reshape((-1, 1))
399 on_plane += line_origins[valid]
401 result = [on_plane, valid]
403 if return_distance:
404 result.append(distance)
405 if return_denom:
406 result.append(projection_dir)
408 return result
411def slice_faces_plane(
412 vertices,
413 faces,
414 plane_normal,
415 plane_origin,
416 uv=None,
417 face_index=None,
418 cached_dots=None,
419):
420 """
421 Slice a mesh (given as a set of faces and vertices) with a plane, returning a
422 new mesh (again as a set of faces and vertices) that is the
423 portion of the original mesh to the positive normal side of the plane.
425 Parameters
426 ---------
427 vertices : (n, 3) float
428 Vertices of source mesh to slice
429 faces : (n, 3) int
430 Faces of source mesh to slice
431 plane_normal : (3,) float
432 Normal vector of plane to intersect with mesh
433 plane_origin : (3,) float
434 Point on plane to intersect with mesh
435 uv : (n, 2) float, optional
436 UV coordinates of source mesh to slice
437 face_index : ((m,) int)
438 Indexes of faces to slice. When no mask is provided, the
439 default is to slice all faces.
440 cached_dots : (n, 3) float
441 If an external function has stored dot
442 products pass them here to avoid recomputing
444 Returns
445 ----------
446 new_vertices : (n, 3) float
447 Vertices of sliced mesh
448 new_faces : (n, 3) int
449 Faces of sliced mesh
450 new_uv : (n, 2) int or None
451 UV coordinates of sliced mesh
452 """
454 if len(vertices) == 0:
455 return vertices, faces, uv
457 have_uv = uv is not None
459 # Construct a mask for the faces to slice.
460 if face_index is not None:
461 faces = faces[face_index]
463 if cached_dots is not None:
464 dots = cached_dots
465 else:
466 # dot product of each vertex with the plane normal indexed by face
467 # so for each face the dot product of each vertex is a row
468 # shape is the same as faces (n,3)
469 dots = np.dot(vertices - plane_origin, plane_normal)
471 # Find vertex orientations w.r.t. faces for all triangles:
472 # -1 -> vertex "inside" plane (positive normal direction)
473 # 0 -> vertex on plane
474 # 1 -> vertex "outside" plane (negative normal direction)
475 signs = np.zeros(len(vertices), dtype=np.int8)
476 signs[dots < -tol.merge] = 1
477 signs[dots > tol.merge] = -1
478 signs = signs[faces]
480 # Find all triangles that intersect this plane
481 # onedge <- indices of all triangles intersecting the plane
482 # inside <- indices of all triangles "inside" the plane (positive normal)
483 signs_sum = signs.sum(axis=1, dtype=np.int8)
484 signs_asum = np.abs(signs).sum(axis=1, dtype=np.int8)
486 # Cases:
487 # (0,0,0), (-1,0,0), (-1,-1,0), (-1,-1,-1) <- inside
488 # (1,0,0), (1,1,0), (1,1,1) <- outside
489 # (1,0,-1), (1,-1,-1), (1,1,-1) <- onedge
490 onedge = np.logical_and(signs_asum >= 2, np.abs(signs_sum) <= 1)
492 inside = signs_sum == -signs_asum
494 # for any faces that lie exactly on-the-plane
495 # we want to only include them if their normal
496 # is backwards from the slicing normal
497 on_plane = signs_asum == 0
498 if on_plane.any():
499 # compute the normals and whether
500 # face is degenerate here
501 check, valid = tm.normals(vertices[faces[on_plane]])
502 # only include faces back from normal
503 dot_check = np.dot(check, plane_normal)
504 # exclude any degenerate faces from the result
505 inside[on_plane] = valid
506 # exclude the degenerate face from our mask
507 on_plane[on_plane] = valid
508 # apply results for this subset
509 inside[on_plane] = dot_check < 0.0
511 # Automatically include all faces that are "inside"
512 new_faces = faces[inside]
514 # Separate faces on the edge into two cases: those which will become
515 # quads (two vertices inside plane) and those which will become triangles
516 # (one vertex inside plane)
517 triangles = vertices[faces]
518 cut_triangles = triangles[onedge]
519 cut_faces_quad = faces[np.logical_and(onedge, signs_sum < 0)]
520 cut_faces_tri = faces[np.logical_and(onedge, signs_sum >= 0)]
521 cut_signs_quad = signs[np.logical_and(onedge, signs_sum < 0)]
522 cut_signs_tri = signs[np.logical_and(onedge, signs_sum >= 0)]
524 # If no faces to cut, the surface is not in contact with this plane.
525 # Thus, return a mesh with only the inside faces
526 if len(cut_faces_quad) + len(cut_faces_tri) == 0:
527 if len(new_faces) == 0:
528 # if no new faces at all return empty arrays
529 empty = (
530 np.zeros((0, 3), dtype=np.float64),
531 np.zeros((0, 3), dtype=np.int64),
532 np.zeros((0, 2), dtype=np.float64) if have_uv else None,
533 )
534 return empty
536 # find the unique indices in the new faces
537 # using an integer-only unique function
538 unique, inverse = grouping.unique_bincount(
539 new_faces.reshape(-1), minlength=len(vertices), return_inverse=True
540 )
542 # use the unique indices for our final vertices and faces
543 final_vert = vertices[unique]
544 final_face = inverse.reshape((-1, 3))
545 final_uv = uv[unique] if have_uv else None
547 return final_vert, final_face, final_uv
549 # Extract the intersections of each triangle's edges with the plane
550 o = cut_triangles # origins
551 d = np.roll(o, -1, axis=1) - o # directions
552 num = (plane_origin - o).dot(plane_normal) # compute num/denom
553 denom = np.dot(d, plane_normal)
554 denom[denom == 0.0] = 1e-12 # prevent division by zero
555 dist = np.divide(num, denom)
556 # intersection points for each segment
557 int_points = dist[:, :, None] * d + o
559 # Initialize the array of new vertices with the current vertices
560 new_vertices = vertices
561 new_quad_vertices = np.zeros((0, 3))
562 new_tri_vertices = np.zeros((0, 3))
564 # Handle the case where a new quad is formed by the intersection
565 # First, extract the intersection points belonging to a new quad
566 quad_int_points = int_points[(signs_sum < 0)[onedge], :, :]
567 num_quads = len(quad_int_points)
568 if num_quads > 0:
569 # Extract the vertex on the outside of the plane, then get the vertices
570 # (in CCW order of the inside vertices)
571 quad_int_inds = np.where(cut_signs_quad == 1)[1]
572 quad_int_verts = cut_faces_quad[
573 np.stack((range(num_quads), range(num_quads)), axis=1),
574 np.stack(((quad_int_inds + 1) % 3, (quad_int_inds + 2) % 3), axis=1),
575 ]
577 # Fill out new quad faces with the intersection points as vertices
578 new_quad_faces = np.append(
579 quad_int_verts,
580 np.arange(len(new_vertices), len(new_vertices) + 2 * num_quads).reshape(
581 num_quads, 2
582 ),
583 axis=1,
584 )
586 # Extract correct intersection points from int_points and order them in
587 # the same way as they were added to faces
588 new_quad_vertices = quad_int_points[
589 np.stack((range(num_quads), range(num_quads)), axis=1),
590 np.stack((((quad_int_inds + 2) % 3).T, quad_int_inds.T), axis=1),
591 :,
592 ].reshape(2 * num_quads, 3)
594 # Add new vertices to existing vertices, triangulate quads, and add the
595 # resulting triangles to the new faces
596 new_vertices = np.append(new_vertices, new_quad_vertices, axis=0)
597 new_tri_faces_from_quads = geometry.triangulate_quads(new_quad_faces)
598 new_faces = np.append(new_faces, new_tri_faces_from_quads, axis=0)
600 # Handle the case where a new triangle is formed by the intersection
601 # First, extract the intersection points belonging to a new triangle
602 tri_int_points = int_points[(signs_sum >= 0)[onedge], :, :]
603 num_tris = len(tri_int_points)
604 if num_tris > 0:
605 # Extract the single vertex for each triangle inside the plane and get the
606 # inside vertices (CCW order)
607 tri_int_inds = np.where(cut_signs_tri == -1)[1]
608 tri_int_verts = cut_faces_tri[range(num_tris), tri_int_inds].reshape(num_tris, 1)
610 # Fill out new triangles with the intersection points as vertices
611 new_tri_faces = np.append(
612 tri_int_verts,
613 np.arange(len(new_vertices), len(new_vertices) + 2 * num_tris).reshape(
614 num_tris, 2
615 ),
616 axis=1,
617 )
619 # Extract correct intersection points and order them in the same way as
620 # the vertices were added to the faces
621 new_tri_vertices = tri_int_points[
622 np.stack((range(num_tris), range(num_tris)), axis=1),
623 np.stack((tri_int_inds.T, ((tri_int_inds + 2) % 3).T), axis=1),
624 :,
625 ].reshape(2 * num_tris, 3)
627 # Append new vertices and new faces
628 new_vertices = np.append(new_vertices, new_tri_vertices, axis=0)
629 new_faces = np.append(new_faces, new_tri_faces, axis=0)
631 # find the unique indices in the new faces
632 # using an integer-only unique function
633 unique, inverse = grouping.unique_bincount(
634 new_faces.reshape(-1), minlength=len(new_vertices), return_inverse=True
635 )
637 # use the unique indexes for our final vertex and faces
638 final_vert = new_vertices[unique]
639 final_face = inverse.reshape((-1, 3))
641 final_uv = None
642 if have_uv:
643 # Generate barycentric coordinates for intersection vertices
644 quad_barycentrics = points_to_barycentric(
645 np.repeat(vertices[cut_faces_quad], 2, axis=0), new_quad_vertices
646 )
647 tri_barycentrics = points_to_barycentric(
648 np.repeat(vertices[cut_faces_tri], 2, axis=0), new_tri_vertices
649 )
650 all_barycentrics = np.concatenate([quad_barycentrics, tri_barycentrics])
652 # Interpolate UVs
653 cut_uv = np.concatenate([uv[cut_faces_quad], uv[cut_faces_tri]])
654 new_uv = np.einsum("ijk,ij->ik", np.repeat(cut_uv, 2, axis=0), all_barycentrics)
655 final_uv = np.concatenate([uv, new_uv])[unique]
657 return final_vert, final_face, final_uv
660def slice_mesh_plane(
661 mesh,
662 plane_normal,
663 plane_origin,
664 face_index=None,
665 cap=False,
666 engine=None,
667 **kwargs,
668):
669 """
670 Slice a mesh with a plane returning a new mesh that is the
671 portion of the original mesh to the positive normal side
672 of the plane.
674 Parameters
675 ---------
676 mesh : Trimesh object
677 Source mesh to slice
678 plane_normal : (3,) float
679 Normal vector of plane to intersect with mesh
680 plane_origin : (3,) float
681 Point on plane to intersect with mesh
682 cap : bool
683 If True, cap the result with a triangulated polygon
684 face_index : ((m,) int)
685 Indexes of mesh.faces to slice. When no mask is provided, the
686 default is to slice all faces.
687 cached_dots : (n, 3) float
688 If an external function has stored dot
689 products pass them here to avoid recomputing
690 engine : None or str
691 Triangulation engine passed to `triangulate_polygon`
692 kwargs : dict
693 Passed to the newly created sliced mesh
695 Returns
696 ----------
697 new_mesh : Trimesh object
698 Sliced mesh
699 """
700 # check input for none
701 if mesh is None:
702 return None
704 # avoid circular import
705 from scipy.spatial import cKDTree
707 from .base import Trimesh
708 from .creation import triangulate_polygon
709 from .path import polygons
710 from .visual import TextureVisuals
712 # check input plane
713 plane_normal = np.asanyarray(plane_normal, dtype=np.float64)
714 plane_origin = np.asanyarray(plane_origin, dtype=np.float64)
716 # check to make sure origins and normals have acceptable shape
717 shape_ok = (
718 (plane_origin.shape == (3,) or util.is_shape(plane_origin, (-1, 3)))
719 and (plane_normal.shape == (3,) or util.is_shape(plane_normal, (-1, 3)))
720 and plane_origin.shape == plane_normal.shape
721 )
722 if not shape_ok:
723 raise ValueError("plane origins and normals must be (n, 3)!")
725 # start with copy of original mesh, faces, and vertices
726 vertices = mesh.vertices.copy()
727 faces = mesh.faces.copy()
729 # We copy the UV coordinates if available
730 has_uv = (
731 hasattr(mesh.visual, "uv") and np.shape(mesh.visual.uv) == (len(mesh.vertices), 2)
732 ) and not cap
733 uv = mesh.visual.uv.copy() if has_uv else None
735 if "process" not in kwargs:
736 kwargs["process"] = False
738 # slice away specified planes
739 for origin, normal in zip(
740 plane_origin.reshape((-1, 3)), plane_normal.reshape((-1, 3))
741 ):
742 # save the new vertices and faces
743 vertices, faces, uv = slice_faces_plane(
744 vertices=vertices,
745 faces=faces,
746 uv=uv,
747 plane_normal=normal,
748 plane_origin=origin,
749 face_index=face_index,
750 )
751 # check if cap arg specified
752 if cap:
753 if face_index:
754 # This hasn't been implemented yet.
755 raise NotImplementedError("face_index and cap can't be used together")
757 # start by deduplicating vertices again
758 unique, inverse = grouping.unique_rows(vertices)
759 vertices = vertices[unique]
760 # will collect additional faces
761 f = inverse[faces]
762 # remove degenerate faces by checking to make sure
763 # that each face has three unique indices
764 f = f[(f[:, :1] != f[:, 1:]).all(axis=1)]
765 # transform to the cap plane
766 to_2D = geometry.plane_transform(origin=origin, normal=-normal)
767 to_3D = np.linalg.inv(to_2D)
769 vertices_2D = tf.transform_points(vertices, to_2D)
770 edges = geometry.faces_to_edges(f)
771 edges.sort(axis=1)
773 on_plane = np.abs(vertices_2D[:, 2]) < 1e-8
774 edges = edges[on_plane[edges].all(axis=1)]
775 edges = edges[edges[:, 0] != edges[:, 1]]
777 unique_edge = grouping.group_rows(edges, require_count=1)
778 if len(unique) < 3:
779 continue
781 tree = cKDTree(vertices)
782 # collect new faces
783 faces = [f]
784 for p in polygons.edges_to_polygons(edges[unique_edge], vertices_2D[:, :2]):
785 # triangulate cap and raise an error if any new vertices were inserted
786 vn, fn = triangulate_polygon(p, engine=engine, force_vertices=True)
787 # collect the original index for the new vertices
788 vn3 = tf.transform_points(util.stack_3D(vn), to_3D)
789 distance, vid = tree.query(vn3)
790 if distance.max() > 1e-8:
791 util.log.debug("triangulate may have inserted vertex!")
792 # triangulation should not have inserted vertices
793 nf = vid[fn]
794 # hmm but it may have returned faces that are now degenerate
795 nf_ok = (nf[:, 1:] != nf[:, :1]).all(axis=1) & (nf[:, 1] != nf[:, 2])
796 faces.append(nf[nf_ok])
798 faces = np.vstack(faces)
800 visual = (
801 TextureVisuals(uv=uv, material=mesh.visual.material.copy()) if has_uv else None
802 )
804 # return the sliced mesh
805 return Trimesh(vertices=vertices, faces=faces, visual=visual, **kwargs)