Coverage for trimesh/geometry.py: 86%
153 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 NDArray
7try:
8 import scipy.sparse
9except BaseException as E:
10 from . import exceptions
12 # raise E again if anyone tries to use sparse
13 scipy = exceptions.ExceptionWrapper(E)
16def plane_transform(origin, normal):
17 """
18 Given the origin and normal of a plane find the transform
19 that will move that plane to be coplanar with the XY plane.
21 Parameters
22 ----------
23 origin : (3,) float
24 Point that lies on the plane
25 normal : (3,) float
26 Vector that points along normal of plane
28 Returns
29 ---------
30 transform: (4,4) float
31 Transformation matrix to move points onto XY plane
32 """
33 transform = align_vectors(normal, [0, 0, 1])
34 if origin is not None:
35 transform[:3, 3] = -np.dot(transform, np.append(origin, 1))[:3]
36 return transform
39def align_vectors(a, b, return_angle=False):
40 """
41 Find the rotation matrix that transforms one 3D vector
42 to another.
44 Parameters
45 ------------
46 a : (3,) float
47 Unit vector
48 b : (3,) float
49 Unit vector
50 return_angle : bool
51 Return the angle between vectors or not
53 Returns
54 -------------
55 matrix : (4, 4) float
56 Homogeneous transform to rotate from `a` to `b`
57 angle : float
58 If `return_angle` angle in radians between `a` and `b`
60 """
61 a = np.array(a, dtype=np.float64)
62 b = np.array(b, dtype=np.float64)
63 if a.shape != (3,) or b.shape != (3,):
64 raise ValueError("vectors must be (3,)!")
66 # find the SVD of the two vectors
67 au = np.linalg.svd(a.reshape((-1, 1)))[0]
68 bu = np.linalg.svd(b.reshape((-1, 1)))[0]
70 if np.linalg.det(au) < 0:
71 au[:, -1] *= -1.0
72 if np.linalg.det(bu) < 0:
73 bu[:, -1] *= -1.0
75 # put rotation into homogeneous transformation
76 matrix = np.eye(4)
77 matrix[:3, :3] = bu.dot(au.T)
79 if return_angle:
80 # projection of a onto b
81 # first row of SVD result is normalized source vector
82 dot = np.dot(au[0], bu[0])
83 # clip to avoid floating point error
84 angle = np.arccos(np.clip(dot, -1.0, 1.0))
85 if dot < -1e-5:
86 angle += np.pi
87 return matrix, angle
89 return matrix
92def faces_to_edges(faces, return_index=False):
93 """
94 Given a list of faces (n,3), return a list of edges (n*3,2)
96 Parameters
97 -----------
98 faces : (n, 3) int
99 Vertex indices representing faces
101 Returns
102 -----------
103 edges : (n*3, 2) int
104 Vertex indices representing edges
105 """
106 faces = np.asanyarray(faces, np.int64)
108 # each face has three edges
109 edges = faces[:, [0, 1, 1, 2, 2, 0]].reshape((-1, 2))
111 if return_index:
112 # edges are in order of faces due to reshape
113 face_index = np.tile(np.arange(len(faces)), (3, 1)).T.reshape(-1)
114 return edges, face_index
115 return edges
118def vector_angle(pairs):
119 """
120 Find the angles between pairs of unit vectors.
122 Parameters
123 ----------
124 pairs : (n, 2, 3) float
125 Unit vector pairs
127 Returns
128 ----------
129 angles : (n,) float
130 Angles between vectors in radians
131 """
132 pairs = np.asanyarray(pairs, dtype=np.float64)
133 if len(pairs) == 0:
134 return np.array([])
135 elif util.is_shape(pairs, (2, 3)):
136 pairs = pairs.reshape((-1, 2, 3))
137 elif not util.is_shape(pairs, (-1, 2, (2, 3))):
138 raise ValueError("pairs must be (n,2,(2|3))!")
140 # do the dot product between vectors
141 dots = util.diagonal_dot(pairs[:, 0], pairs[:, 1])
142 # clip for floating point error
143 dots = np.clip(dots, -1.0, 1.0)
144 # do cos and remove arbitrary sign
145 angles = np.abs(np.arccos(dots))
147 return angles
150def triangulate_quads(quads, dtype=np.int64, use_fan: bool = True) -> NDArray:
151 """
152 Given an array of quad faces return them as triangle faces,
153 also handles pure triangles and mixed triangles and quads.
155 Parameters
156 -----------
157 quads: (n, 4) int
158 Vertex indices of quad faces.
159 dtype
160 Data type requested for the return
161 use_fan
162 Triangulate holes larger than quads with fans,
163 which may be wrong if the holes are non-convex
165 Returns
166 -----------
167 faces : (m, 3) int
168 Vertex indices of triangular faces.
169 """
171 if len(quads) == 0:
172 return np.zeros(0, dtype=dtype)
174 try:
175 # this will fail in newer versions of numpy
176 # if there are mixed quads and tris
177 quads = np.array(quads, dtype=dtype)
179 if len(quads.shape) == 2 and quads.shape[1] == 3:
180 # if they are just triangles return immediately
181 return quads.astype(dtype)
183 if len(quads.shape) == 2 and quads.shape[1] == 4:
184 # if they are just quads stack and return
185 return np.vstack((quads[:, [0, 1, 2]], quads[:, [2, 3, 0]])).astype(dtype)
186 except ValueError:
187 # new numpy raises an error for sequences
188 pass
190 # we made it here so we have mixed tris/quads/polygons
191 # do one pass to get the lengths
192 lengths = np.array([len(i) for i in quads], dtype=np.int64)
194 # get triangles and quads as clean constant-row numpy arrays
195 tri = np.array([quads[i] for i in np.nonzero(lengths == 3)[0]], dtype=np.int64)
196 quad = np.array([quads[i] for i in np.nonzero(lengths == 4)[0]], dtype=np.int64)
198 if use_fan:
199 # get arbitrary polygons as a ragged sequence
200 poly = [quads[i] for i in np.nonzero(lengths > 4)[0]]
201 else:
202 # skip any hole larger than a quad
203 poly = []
205 if len(quad) == 0 and len(poly) == 0:
206 # only triangles, return triangles
207 return tri.astype(dtype)
208 if len(poly) > 0:
209 # use numpy slicing to triangulate ragged-sequence polygons
210 poly = util.triangle_fans_to_faces(poly)
211 if len(quad) > 0:
212 # trivially tessellate quads
213 quad = np.vstack((quad[:, [0, 1, 2]], quad[:, [2, 3, 0]]))
214 # stack triangles from all three cases
215 return util.vstack_empty([tri, quad, poly]).astype(dtype)
218def vertex_face_indices(vertex_count, faces, faces_sparse):
219 """
220 Find vertex face indices from the faces array of vertices
222 Parameters
223 -----------
224 vertex_count : int
225 The number of vertices faces refer to
226 faces : (n, 3) int
227 List of vertex indices
228 faces_sparse : scipy.sparse.COO
229 Sparse matrix
231 Returns
232 -----------
233 vertex_faces : (vertex_count, ) int
234 Face indices for every vertex
235 Array padded with -1 in each row for all vertices with fewer
236 face indices than the max number of face indices.
237 """
238 # Create 2D array with row for each vertex and
239 # length of max number of faces for a vertex
240 try:
241 counts = np.bincount(faces.flatten(), minlength=vertex_count)
242 except TypeError:
243 # casting failed on 32 bit Windows
244 log.warning("casting failed, falling back!")
245 # fall back to np.unique (usually ~35x slower than bincount)
246 counts = np.unique(faces.flatten(), return_counts=True)[1]
247 assert len(counts) == vertex_count
248 assert faces.max() < vertex_count
250 # start cumulative sum at zero and clip off the last value
251 starts = np.append(0, np.cumsum(counts)[:-1])
252 # pack incrementing array into final shape
253 pack = np.arange(counts.max()) + starts[:, None]
254 # pad each row with -1 to pad to the max length
255 padded = -(pack >= (starts + counts)[:, None]).astype(np.int64)
257 try:
258 # do most of the work with a sparse dot product
259 identity = scipy.sparse.identity(len(faces), dtype=int)
260 sorted_faces = faces_sparse.dot(identity).nonzero()[1]
261 # this will fail if any face was degenerate
262 # TODO
263 # figure out how to filter out degenerate faces from sparse
264 # result if sorted_faces.size != faces.size
265 padded[padded == 0] = sorted_faces
266 except BaseException:
267 # fall back to a slow loop
268 log.warning(
269 "vertex_faces falling back to slow loop! "
270 + "mesh probably has degenerate faces",
271 exc_info=True,
272 )
273 sort = np.zeros(faces.size, dtype=np.int64)
274 flat = faces.flatten()
275 for v in range(vertex_count):
276 # assign the data in order
277 sort[starts[v] : starts[v] + counts[v]] = (np.where(flat == v)[0] // 3)[::-1]
278 padded[padded == 0] = sort
279 return padded
282def mean_vertex_normals(vertex_count, faces, face_normals, sparse=None, **kwargs):
283 """
284 Find vertex normals from the mean of the faces that contain
285 that vertex.
287 Parameters
288 -----------
289 vertex_count : int
290 The number of vertices faces refer to
291 faces : (n, 3) int
292 List of vertex indices
293 face_normals : (n, 3) float
294 Normal vector for each face
296 Returns
297 -----------
298 vertex_normals : (vertex_count, 3) float
299 Normals for every vertex
300 Vertices unreferenced by faces will be zero.
301 """
303 def summed_sparse():
304 # use a sparse matrix of which face contains each vertex to
305 # figure out the summed normal at each vertex
306 # allow cached sparse matrix to be passed
307 if sparse is None:
308 matrix = index_sparse(vertex_count, faces)
309 else:
310 matrix = sparse
311 summed = matrix.dot(face_normals)
312 return summed
314 def summed_loop():
315 # loop through every face, in tests was ~50x slower than
316 # doing this with a sparse matrix
317 summed = np.zeros((vertex_count, 3))
318 for face, normal in zip(faces, face_normals):
319 summed[face] += normal
320 return summed
322 try:
323 summed = summed_sparse()
324 except BaseException:
325 log.warning("unable to use sparse matrix, falling back!", exc_info=True)
326 summed = summed_loop()
328 # invalid normals will be returned as zero
329 vertex_normals = util.unitize(summed)
331 return vertex_normals
334def weighted_vertex_normals(
335 vertex_count, faces, face_normals, face_angles, use_loop=False
336):
337 """
338 Compute vertex normals from the faces that contain that vertex.
339 The contribution of a face's normal to a vertex normal is the
340 ratio of the corner-angle in which the vertex is, with respect
341 to the sum of all corner-angles surrounding the vertex.
343 Grit Thuerrner & Charles A. Wuethrich (1998)
344 Computing Vertex Normals from Polygonal Facets,
345 Journal of Graphics Tools, 3:1, 43-46
347 Parameters
348 -----------
349 vertex_count : int
350 The number of vertices faces refer to
351 faces : (n, 3) int
352 List of vertex indices
353 face_normals : (n, 3) float
354 Normal vector for each face
355 face_angles : (n, 3) float
356 Angles at each vertex in the face
358 Returns
359 -----------
360 vertex_normals : (vertex_count, 3) float
361 Normals for every vertex
362 Vertices unreferenced by faces will be zero.
363 """
365 def summed_sparse():
366 # use a sparse matrix of which face contains each vertex to
367 # figure out the summed normal at each vertex
368 # allow cached sparse matrix to be passed
369 # fill the matrix with vertex-corner angles as weights
370 matrix = index_sparse(vertex_count, faces, data=face_angles.ravel())
371 return matrix.dot(face_normals)
373 def summed_loop():
374 summed = np.zeros((vertex_count, 3), np.float64)
375 for vertex_idx in np.arange(vertex_count):
376 # loop over all vertices
377 # compute normal contributions from surrounding faces
378 # obviously slower than with the sparse matrix
379 face_idxs, inface_idxs = np.where(faces == vertex_idx)
380 surrounding_angles = face_angles[face_idxs, inface_idxs]
381 summed[vertex_idx] = np.dot(
382 surrounding_angles / surrounding_angles.sum(), face_normals[face_idxs]
383 )
385 return summed
387 # normals should be unit vectors
388 face_ok = (face_normals**2).sum(axis=1) > 0.5
389 # don't consider faces with invalid normals
390 faces = faces[face_ok]
391 face_normals = face_normals[face_ok]
392 face_angles = face_angles[face_ok]
394 if not use_loop:
395 try:
396 return util.unitize(summed_sparse())
397 except BaseException:
398 log.warning("unable to use sparse matrix, falling back!", exc_info=True)
399 # we either crashed or were asked to loop
400 return util.unitize(summed_loop())
403def index_sparse(columns, indices, data=None, dtype=None):
404 """
405 Return a sparse matrix for which vertices are contained in which faces.
406 A data vector can be passed which is then used instead of booleans
408 Parameters
409 ------------
410 columns : int
411 Number of columns, usually number of vertices
412 indices : (m, d) int
413 Usually mesh.faces
415 Returns
416 ---------
417 sparse: scipy.sparse.coo_matrix of shape (columns, len(faces))
418 dtype is boolean
420 Examples
421 ----------
422 In [1]: sparse = faces_sparse(len(mesh.vertices), mesh.faces)
424 In [2]: sparse.shape
425 Out[2]: (12, 20)
427 In [3]: mesh.faces.shape
428 Out[3]: (20, 3)
430 In [4]: mesh.vertices.shape
431 Out[4]: (12, 3)
433 In [5]: dense = sparse.toarray().astype(int)
435 In [6]: dense
436 Out[6]:
437 array([[1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
438 [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
439 [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
440 [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0],
441 [0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
442 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0],
443 [0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
444 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1],
445 [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0],
446 [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0],
447 [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1],
448 [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1]])
450 In [7]: dense.sum(axis=0)
451 Out[7]: array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3])
452 """
453 indices = np.asanyarray(indices)
454 columns = int(columns)
456 # flattened list
457 row = indices.reshape(-1)
458 col = np.tile(
459 np.arange(len(indices)).reshape((-1, 1)), (1, indices.shape[1])
460 ).reshape(-1)
462 shape = (columns, len(indices))
463 if data is None:
464 data = np.ones(len(col), dtype=bool)
465 elif len(data) != len(col):
466 raise ValueError("data incorrect length")
468 if dtype is not None:
469 data = data.astype(dtype)
471 # assemble into sparse matrix
472 return scipy.sparse.coo_matrix((data, (row, col)), shape=shape, dtype=data.dtype)