Coverage for trimesh/triangles.py: 86%
262 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2triangles.py
3-------------
5Functions for dealing with triangle soups in (n, 3, 3) float form.
6"""
8from dataclasses import dataclass
9from logging import getLogger
11import numpy as np
13from . import util
14from .constants import tol
15from .points import point_plane_distance
16from .typed import NDArray, float64
17from .util import diagonal_dot, unitize
19log = getLogger(__name__)
22def cross(triangles: NDArray) -> NDArray:
23 """
24 Returns the cross product of two edges from input triangles
26 Parameters
27 --------------
28 triangles: (n, 3, 3) float
29 Vertices of triangles
31 Returns
32 --------------
33 crosses : (n, 3) float
34 Cross product of two edge vectors
35 """
36 vectors = triangles[:, 1:, :] - triangles[:, :2, :]
37 if triangles.shape[2] == 3:
38 return np.cross(vectors[:, 0], vectors[:, 1])
39 elif triangles.shape[2] == 2:
40 a = vectors[:, 0]
41 b = vectors[:, 1]
42 # numpy 2.0 deprecated 2D cross productes
43 return a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0]
45 raise ValueError(triangles.shape)
48def area(triangles=None, crosses=None):
49 """
50 Calculates the sum area of input triangles
52 Parameters
53 ----------
54 triangles : (n, 3, 3) float
55 Vertices of triangles
56 crosses : (n, 3) float or None
57 As a speedup don't re- compute cross products
58 sum : bool
59 Return summed area or individual triangle area
61 Returns
62 ----------
63 area : (n,) float or float
64 Individual or summed area depending on `sum` argument
65 """
66 if crosses is None:
67 crosses = cross(np.asanyarray(triangles, dtype=np.float64))
68 if len(crosses.shape) == 1:
69 # support 2D triangles
70 return np.abs(crosses) / 2.0
71 return np.sqrt((crosses**2).sum(axis=1)) / 2.0
74def normals(triangles=None, crosses=None):
75 """
76 Calculates the normals of input triangles
78 Parameters
79 ------------
80 triangles : (n, 3, 3) float
81 Vertex positions
82 crosses : (n, 3) float
83 Cross products of edge vectors
85 Returns
86 ------------
87 normals : (m, 3) float
88 Normal vectors
89 valid : (n,) bool
90 Was the face nonzero area or not
91 """
92 if triangles is not None:
93 triangles = np.asanyarray(triangles, dtype=np.float64)
94 if triangles.shape[-1] == 2:
95 # 2D triangles return unit normal along Z
96 unit = np.tile([0.0, 0.0, 1.0], (triangles.shape[0], 1))
97 valid = np.ones(len(triangles), dtype=bool)
98 return unit, valid
99 if crosses is None:
100 crosses = cross(triangles)
101 # unitize the cross product vectors
102 unit, valid = unitize(crosses, check_valid=True)
103 return unit, valid
106def angles(triangles):
107 """
108 Calculates the angles of input triangles.
110 Parameters
111 ------------
112 triangles : (n, 3, 3) float
113 Vertex positions
115 Returns
116 ------------
117 angles : (n, 3) float
118 Angles at vertex positions in radians
119 Degenerate angles will be returned as zero
120 """
121 # don't copy triangles
122 triangles = np.asanyarray(triangles, dtype=np.float64)
124 # get a unit vector for each edge of the triangle
125 u = unitize(triangles[:, 1] - triangles[:, 0])
126 v = unitize(triangles[:, 2] - triangles[:, 0])
127 w = unitize(triangles[:, 2] - triangles[:, 1])
129 # run the cosine and per-row dot product
130 result = np.zeros((len(triangles), 3), dtype=np.float64)
131 # clip to make sure we don't float error past 1.0
132 result[:, 0] = np.arccos(np.clip(diagonal_dot(u, v), -1, 1))
133 result[:, 1] = np.arccos(np.clip(diagonal_dot(-u, w), -1, 1))
134 # the third angle is just the remaining
135 result[:, 2] = np.pi - result[:, 0] - result[:, 1]
137 # a triangle with any zero angles is degenerate
138 # so set all of the angles to zero in that case
139 result[(result < tol.merge).any(axis=1), :] = 0.0
141 return result
144def all_coplanar(triangles):
145 """
146 Check to see if a list of triangles are all coplanar
148 Parameters
149 ----------------
150 triangles: (n, 3, 3) float
151 Vertices of triangles
153 Returns
154 ---------------
155 all_coplanar : bool
156 True if all triangles are coplanar
157 """
158 triangles = np.asanyarray(triangles, dtype=np.float64)
159 if not util.is_shape(triangles, (-1, 3, 3)):
160 raise ValueError("Triangles must be (n, 3, 3)!")
162 test_normal = normals(triangles)[0]
163 test_vertex = triangles[0][0]
164 distances = point_plane_distance(
165 points=triangles[1:].reshape((-1, 3)),
166 plane_normal=test_normal,
167 plane_origin=test_vertex,
168 )
169 all_coplanar = np.all(np.abs(distances) < tol.zero)
170 return all_coplanar
173def any_coplanar(triangles):
174 """
175 For a list of triangles if the FIRST triangle is coplanar
176 with ANY of the following triangles, return True.
177 Otherwise, return False.
178 """
179 triangles = np.asanyarray(triangles, dtype=np.float64)
180 if not util.is_shape(triangles, (-1, 3, 3)):
181 raise ValueError("Triangles must be (n, 3, 3)!")
183 test_normal = normals(triangles)[0]
184 test_vertex = triangles[0][0]
185 distances = point_plane_distance(
186 points=triangles[1:].reshape((-1, 3)),
187 plane_normal=test_normal,
188 plane_origin=test_vertex,
189 )
190 any_coplanar = np.any(np.all(np.abs(distances.reshape((-1, 3)) < tol.zero), axis=1))
191 return any_coplanar
194@dataclass
195class MassProperties:
196 # the density value these mass properties were calculated with
197 # this alters `mass` and `inertia`
198 density: float
200 # the volume multiplied by the density
201 mass: float
203 # the volume produced
204 volume: float
206 # the (3,) center of mass
207 center_mass: NDArray[float64]
209 # the (3, 3) inertia tensor
210 inertia: NDArray[float64] | None = None
212 def __getitem__(self, item: str):
213 # add for backwards compatibility
214 return getattr(self, item)
217def mass_properties(
218 triangles, crosses=None, density=None, center_mass=None, skip_inertia=False
219) -> MassProperties:
220 """
221 Calculate the mass properties of a group of triangles.
223 Implemented from:
224 http://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf
226 Parameters
227 ----------
228 triangles : (n, 3, 3) float
229 Triangle vertices in space
230 crosses : (n,) float
231 Optional cross products of triangles
232 density : float
233 Optional override for density
234 center_mass : (3,) float
235 Optional override for center mass
236 skip_inertia : bool
237 if True will not return moments matrix
239 Returns
240 ---------
241 info : dict
242 Mass properties
243 """
244 triangles = np.asanyarray(triangles, dtype=np.float64)
245 if not util.is_shape(triangles, (-1, 3, 3)):
246 raise ValueError("Triangles must be (n, 3, 3)!")
248 if crosses is None:
249 crosses = cross(triangles)
250 if density is None:
251 density = 1.0
253 # these are the subexpressions of the integral
254 # this is equvilant but 7x faster than triangles.sum(axis=1)
255 f1 = triangles[:, 0, :] + triangles[:, 1, :] + triangles[:, 2, :]
257 # for the the first vertex of every triangle:
258 # triangles[:,0,:] will give rows like [[x0, y0, z0], ...]
260 # for the x coordinates of every triangle
261 # triangles[:,:,0] will give rows like [[x0, x1, x2], ...]
262 f2 = (
263 triangles[:, 0, :] ** 2
264 + triangles[:, 1, :] ** 2
265 + triangles[:, 0, :] * triangles[:, 1, :]
266 + triangles[:, 2, :] * f1
267 )
268 f3 = (
269 (triangles[:, 0, :] ** 3)
270 + (triangles[:, 0, :] ** 2) * (triangles[:, 1, :])
271 + (triangles[:, 0, :]) * (triangles[:, 1, :] ** 2)
272 + (triangles[:, 1, :] ** 3)
273 + (triangles[:, 2, :] * f2)
274 )
275 g0 = f2 + (triangles[:, 0, :] + f1) * triangles[:, 0, :]
276 g1 = f2 + (triangles[:, 1, :] + f1) * triangles[:, 1, :]
277 g2 = f2 + (triangles[:, 2, :] + f1) * triangles[:, 2, :]
278 integral = np.zeros((10, len(f1)))
279 integral[0] = crosses[:, 0] * f1[:, 0]
280 integral[1:4] = (crosses * f2).T
281 integral[4:7] = (crosses * f3).T
282 for i in range(3):
283 triangle_i = np.mod(i + 1, 3)
284 integral[i + 7] = crosses[:, i] * (
285 (triangles[:, 0, triangle_i] * g0[:, i])
286 + (triangles[:, 1, triangle_i] * g1[:, i])
287 + (triangles[:, 2, triangle_i] * g2[:, i])
288 )
290 integrated = integral.sum(axis=1) / np.array(
291 [6, 24, 24, 24, 60, 60, 60, 120, 120, 120], dtype=np.float64
292 )
294 volume = integrated[0]
296 # we allow center of mass to be overridden
297 if center_mass is None:
298 if np.abs(volume) < tol.zero:
299 # if there is no volume set center of mass to the centroid
300 log.debug("volume is negative center of mass is almost certain to be wrong!")
301 # otherwise get it from the integration
302 center_mass = integrated[1:4] / volume
304 result = MassProperties(
305 density=density,
306 mass=density * volume,
307 volume=volume,
308 center_mass=center_mass,
309 )
311 if skip_inertia:
312 return result
314 inertia = np.zeros((3, 3))
315 inertia[0, 0] = (
316 integrated[5] + integrated[6] - (volume * (center_mass[[1, 2]] ** 2).sum())
317 )
318 inertia[1, 1] = (
319 integrated[4] + integrated[6] - (volume * (center_mass[[0, 2]] ** 2).sum())
320 )
321 inertia[2, 2] = (
322 integrated[4] + integrated[5] - (volume * (center_mass[[0, 1]] ** 2).sum())
323 )
324 inertia[0, 1] = -(integrated[7] - (volume * np.prod(center_mass[[0, 1]])))
325 inertia[1, 2] = -(integrated[8] - (volume * np.prod(center_mass[[1, 2]])))
326 inertia[0, 2] = -(integrated[9] - (volume * np.prod(center_mass[[0, 2]])))
327 inertia[2, 0] = inertia[0, 2]
328 inertia[2, 1] = inertia[1, 2]
329 inertia[1, 0] = inertia[0, 1]
330 result.inertia = inertia * density
332 return result
335def windings_aligned(triangles, normals_compare):
336 """
337 Given a list of triangles and a list of normals determine if the
338 two are aligned
340 Parameters
341 ----------
342 triangles : (n, 3, 3) float
343 Vertex locations in space
344 normals_compare : (n, 3) float
345 List of normals to compare
347 Returns
348 ----------
349 aligned : (n,) bool
350 Are normals aligned with triangles
351 """
352 triangles = np.asanyarray(triangles, dtype=np.float64)
353 if not util.is_shape(triangles, (-1, 3, 3), allow_zeros=True):
354 raise ValueError(f"triangles must have shape (n, 3, 3), got {triangles.shape!s}")
355 normals_compare = np.asanyarray(normals_compare, dtype=np.float64)
357 calculated, valid = normals(triangles)
358 if normals_compare.shape == (3,):
359 # single comparison vector case
360 difference = np.dot(calculated, normals_compare)
361 else:
362 # multiple comparison case
363 difference = diagonal_dot(calculated, normals_compare[valid])
365 aligned = np.zeros(len(triangles), dtype=bool)
366 aligned[valid] = difference > 0.0
368 return aligned
371def bounds_tree(triangles):
372 """
373 Given a list of triangles, create an r-tree for broad- phase
374 collision detection
376 Parameters
377 ---------
378 triangles : (n, 3, 3) float
379 Triangles in space
381 Returns
382 ---------
383 tree : rtree.Rtree
384 One node per triangle
385 """
386 triangles = np.asanyarray(triangles, dtype=np.float64)
387 if not util.is_shape(triangles, (-1, 3, (2, 3))):
388 raise ValueError("Triangles must be (n, 3, 3)!")
390 # the (n,6) interleaved bounding box for every triangle
391 triangle_bounds = np.column_stack((triangles.min(axis=1), triangles.max(axis=1)))
392 tree = util.bounds_tree(triangle_bounds)
393 return tree
396def nondegenerate(triangles, areas=None, height=None):
397 """
398 Find all triangles which have an oriented bounding box
399 where both of the two sides is larger than a specified height.
401 Degenerate triangles can be when:
402 1) Two of the three vertices are colocated
403 2) All three vertices are unique but colinear
406 Parameters
407 ----------
408 triangles : (n, 3, 3) float
409 Triangles in space
410 height : float
411 Minimum edge length of a triangle to keep
413 Returns
414 ----------
415 nondegenerate : (n,) bool
416 True if a triangle meets required minimum height
417 """
418 triangles = np.asanyarray(triangles, dtype=np.float64)
419 if not util.is_shape(triangles, (-1, 3, 3)):
420 raise ValueError("Triangles must be (n, 3, 3)!")
422 if height is None:
423 height = tol.merge
425 # if both edges of the triangles OBB are longer than tol.merge
426 # we declare them to be nondegenerate
427 ok = (extents(triangles=triangles, areas=areas) > height).all(axis=1)
429 return ok
432def extents(triangles, areas=None):
433 """
434 Return the 2D bounding box size of each triangle.
436 Parameters
437 ----------
438 triangles : (n, 3, 3) float
439 Triangles in space
440 areas : (n,) float
441 Optional area of input triangles
443 Returns
444 ----------
445 box : (n, 2) float
446 The size of each triangle's 2D oriented bounding box
447 """
448 triangles = np.asanyarray(triangles, dtype=np.float64)
449 if not util.is_shape(triangles, (-1, 3, 3)):
450 raise ValueError("Triangles must be (n, 3, 3)!")
452 if areas is None:
453 areas = area(triangles=triangles)
455 # the edge vectors which define the triangle
456 a = triangles[:, 1] - triangles[:, 0]
457 b = triangles[:, 2] - triangles[:, 0]
459 # length of the edge vectors
460 length_a = (a**2).sum(axis=1) ** 0.5
461 length_b = (b**2).sum(axis=1) ** 0.5
463 # which edges are acceptable length
464 nonzero_a = length_a > tol.merge
465 nonzero_b = length_b > tol.merge
467 # find the two heights of the triangle
468 # essentially this is the side length of an
469 # oriented bounding box, per triangle
470 box = np.zeros((len(triangles), 2), dtype=np.float64)
471 box[:, 0][nonzero_a] = (areas[nonzero_a] * 2) / length_a[nonzero_a]
472 box[:, 1][nonzero_b] = (areas[nonzero_b] * 2) / length_b[nonzero_b]
474 return box
477def barycentric_to_points(triangles, barycentric):
478 """
479 Convert a list of barycentric coordinates on a list of triangles
480 to cartesian points.
482 Parameters
483 ------------
484 triangles : (n, 3, 3) float
485 Triangles in space
486 barycentric : (n, 2) float
487 Barycentric coordinates
489 Returns
490 -----------
491 points : (m, 3) float
492 Points in space
493 """
494 barycentric = np.array(barycentric, dtype=np.float64)
495 triangles = np.asanyarray(triangles, dtype=np.float64)
497 # normalize in-place
498 barycentric /= barycentric.sum(axis=1).reshape((-1, 1))
499 points = (triangles * barycentric.reshape((-1, 3, 1))).sum(axis=1)
501 return points
504def points_to_barycentric(triangles, points, method="cramer"):
505 """
506 Find the barycentric coordinates of points relative to triangles.
508 The Cramer's rule solution implements:
509 http://blackpawn.com/texts/pointinpoly
511 The cross product solution implements:
512 https://www.cs.ubc.ca/~heidrich/Papers/JGT.05.pdf
515 Parameters
516 -----------
517 triangles : (n, 3, 2 | 3) float
518 Triangles vertices in space
519 points : (n, 2 | 3) float
520 Point in space associated with a triangle
521 method : str
522 Which method to compute the barycentric coordinates with:
523 - 'cross': uses a method using cross products, roughly 2x slower but
524 different numerical robustness properties
525 - anything else: uses a cramer's rule solution
527 Returns
528 -----------
529 barycentric : (n, 3) float
530 Barycentric coordinates of each point
531 """
533 def method_cross():
534 n = np.cross(edge_vectors[:, 0], edge_vectors[:, 1])
535 denominator = diagonal_dot(n, n)
537 barycentric = np.zeros((len(triangles), 3), dtype=np.float64)
538 barycentric[:, 2] = diagonal_dot(np.cross(edge_vectors[:, 0], w), n) / denominator
539 barycentric[:, 1] = diagonal_dot(np.cross(w, edge_vectors[:, 1]), n) / denominator
540 barycentric[:, 0] = 1 - barycentric[:, 1] - barycentric[:, 2]
541 return barycentric
543 def method_cramer():
544 dot00 = diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 0])
545 dot01 = diagonal_dot(edge_vectors[:, 0], edge_vectors[:, 1])
546 dot02 = diagonal_dot(edge_vectors[:, 0], w)
547 dot11 = diagonal_dot(edge_vectors[:, 1], edge_vectors[:, 1])
548 dot12 = diagonal_dot(edge_vectors[:, 1], w)
550 inverse_denominator = 1.0 / (dot00 * dot11 - dot01 * dot01)
552 barycentric = np.zeros((len(triangles), 3), dtype=np.float64)
553 barycentric[:, 2] = (dot00 * dot12 - dot01 * dot02) * inverse_denominator
554 barycentric[:, 1] = (dot11 * dot02 - dot01 * dot12) * inverse_denominator
555 barycentric[:, 0] = 1 - barycentric[:, 1] - barycentric[:, 2]
556 return barycentric
558 # establish that input triangles and points are sane
559 triangles = np.asanyarray(triangles, dtype=np.float64)
560 points = np.asanyarray(points, dtype=np.float64)
562 # triangles should be (n, 3, dimension)
563 if len(triangles.shape) != 3:
564 raise ValueError("triangles shape incorrect")
566 # this should work for 2D and 3D triangles
567 dim = triangles.shape[2]
568 if (
569 len(points.shape) != 2
570 or points.shape[1] != dim
571 or points.shape[0] != triangles.shape[0]
572 ):
573 raise ValueError("triangles and points must correspond")
575 edge_vectors = triangles[:, 1:] - triangles[:, :1]
576 w = points - triangles[:, 0].reshape((-1, dim))
578 if method == "cross":
579 return method_cross()
580 return method_cramer()
583def closest_point(triangles, points):
584 """
585 Return the closest point on the surface of each triangle for a
586 list of corresponding points.
588 Implements the method from "Real Time Collision Detection" and
589 use the same variable names as "ClosestPtPointTriangle" to avoid
590 being any more confusing.
593 Parameters
594 ----------
595 triangles : (n, 3, 3) float
596 Triangle vertices in space
597 points : (n, 3) float
598 Points in space
600 Returns
601 ----------
602 closest : (n, 3) float
603 Point on each triangle closest to each point
604 """
606 # check input triangles and points
607 triangles = np.asanyarray(triangles, dtype=np.float64)
608 points = np.asanyarray(points, dtype=np.float64)
609 if not util.is_shape(triangles, (-1, 3, 3)):
610 raise ValueError("triangles shape incorrect")
611 if not util.is_shape(points, (len(triangles), 3)):
612 raise ValueError("need same number of triangles and points!")
614 # store the location of the closest point
615 result = np.zeros_like(points)
616 # which points still need to be handled
617 remain = np.ones(len(points), dtype=bool)
619 # if we dot product this against a (n, 3)
620 # it is equivalent but faster than array.sum(axis=1)
621 ones = [1.0, 1.0, 1.0]
623 # get the three points of each triangle
624 # use the same notation as RTCD to avoid confusion
625 a = triangles[:, 0, :]
626 b = triangles[:, 1, :]
627 c = triangles[:, 2, :]
629 # check if P is in vertex region outside A
630 ab = b - a
631 ac = c - a
632 ap = points - a
633 # this is a faster equivalent of:
634 # diagonal_dot(ab, ap)
635 d1 = np.dot(ab * ap, ones)
636 d2 = np.dot(ac * ap, ones)
638 # is the point at A
639 is_a = np.logical_and(d1 < tol.zero, d2 < tol.zero)
640 if any(is_a):
641 result[is_a] = a[is_a]
642 remain[is_a] = False
644 # check if P in vertex region outside B
645 bp = points - b
646 d3 = np.dot(ab * bp, ones)
647 d4 = np.dot(ac * bp, ones)
649 # do the logic check
650 is_b = (d3 > -tol.zero) & (d4 <= d3) & remain
651 if any(is_b):
652 result[is_b] = b[is_b]
653 remain[is_b] = False
655 # check if P in edge region of AB, if so return projection of P onto A
656 vc = (d1 * d4) - (d3 * d2)
657 is_ab = (vc < tol.zero) & (d1 > -tol.zero) & (d3 < tol.zero) & remain
658 if any(is_ab):
659 v = (d1[is_ab] / (d1[is_ab] - d3[is_ab])).reshape((-1, 1))
660 result[is_ab] = a[is_ab] + (v * ab[is_ab])
661 remain[is_ab] = False
663 # check if P in vertex region outside C
664 cp = points - c
665 d5 = np.dot(ab * cp, ones)
666 d6 = np.dot(ac * cp, ones)
667 is_c = (d6 > -tol.zero) & (d5 <= d6) & remain
668 if any(is_c):
669 result[is_c] = c[is_c]
670 remain[is_c] = False
672 # check if P in edge region of AC, if so return projection of P onto AC
673 vb = (d5 * d2) - (d1 * d6)
674 is_ac = (vb < tol.zero) & (d2 > -tol.zero) & (d6 < tol.zero) & remain
675 if any(is_ac):
676 w = (d2[is_ac] / (d2[is_ac] - d6[is_ac])).reshape((-1, 1))
677 result[is_ac] = a[is_ac] + w * ac[is_ac]
678 remain[is_ac] = False
680 # check if P in edge region of BC, if so return projection of P onto BC
681 va = (d3 * d6) - (d5 * d4)
682 is_bc = (va < tol.zero) & ((d4 - d3) > -tol.zero) & ((d5 - d6) > -tol.zero) & remain
683 if any(is_bc):
684 d43 = d4[is_bc] - d3[is_bc]
685 w = (d43 / (d43 + (d5[is_bc] - d6[is_bc]))).reshape((-1, 1))
686 result[is_bc] = b[is_bc] + w * (c[is_bc] - b[is_bc])
687 remain[is_bc] = False
689 # any remaining points must be inside face region
690 if any(remain):
691 # point is inside face region
692 denom = 1.0 / (va[remain] + vb[remain] + vc[remain])
693 v = (vb[remain] * denom).reshape((-1, 1))
694 w = (vc[remain] * denom).reshape((-1, 1))
695 # compute Q through its barycentric coordinates
696 result[remain] = a[remain] + (ab[remain] * v) + (ac[remain] * w)
698 return result
701def to_kwargs(triangles):
702 """
703 Convert a list of triangles to the kwargs for the Trimesh
704 constructor.
706 Parameters
707 ---------
708 triangles : (n, 3, 3) float
709 Triangles in space
711 Returns
712 ---------
713 kwargs : dict
714 Keyword arguments for the trimesh.Trimesh constructor
715 Includes keys 'vertices' and 'faces'
717 Examples
718 ---------
719 >>> mesh = trimesh.Trimesh(**trimesh.triangles.to_kwargs(triangles))
720 """
721 triangles = np.asanyarray(triangles, dtype=np.float64)
722 if not util.is_shape(triangles, (-1, 3, 3)):
723 raise ValueError("Triangles must be (n, 3, 3)!")
725 vertices = triangles.reshape((-1, 3))
726 faces = np.arange(len(vertices)).reshape((-1, 3))
727 kwargs = {"vertices": vertices, "faces": faces}
729 return kwargs