Coverage for trimesh/geometry.py: 86%
151 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-10-02 20:54 +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 column 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 return matrix, angle
87 return matrix
90def faces_to_edges(faces, return_index=False):
91 """
92 Given a list of faces (n,3), return a list of edges (n*3,2)
94 Parameters
95 -----------
96 faces : (n, 3) int
97 Vertex indices representing faces
99 Returns
100 -----------
101 edges : (n*3, 2) int
102 Vertex indices representing edges
103 """
104 faces = np.asanyarray(faces, np.int64)
106 # each face has three edges
107 edges = faces[:, [0, 1, 1, 2, 2, 0]].reshape((-1, 2))
109 if return_index:
110 # edges are in order of faces due to reshape
111 face_index = np.tile(np.arange(len(faces)), (3, 1)).T.reshape(-1)
112 return edges, face_index
113 return edges
116def vector_angle(pairs):
117 """
118 Find the angles between pairs of unit vectors.
120 Parameters
121 ----------
122 pairs : (n, 2, 3) float
123 Unit vector pairs
125 Returns
126 ----------
127 angles : (n,) float
128 Angles between vectors in radians
129 """
130 pairs = np.asanyarray(pairs, dtype=np.float64)
131 if len(pairs) == 0:
132 return np.array([])
133 elif util.is_shape(pairs, (2, 3)):
134 pairs = pairs.reshape((-1, 2, 3))
135 elif not util.is_shape(pairs, (-1, 2, (2, 3))):
136 raise ValueError("pairs must be (n,2,(2|3))!")
138 # do the dot product between vectors
139 dots = util.diagonal_dot(pairs[:, 0], pairs[:, 1])
140 # clip for floating point error
141 dots = np.clip(dots, -1.0, 1.0)
142 # do cos and remove arbitrary sign
143 angles = np.abs(np.arccos(dots))
145 return angles
148def triangulate_quads(quads, dtype=np.int64, use_fan: bool = True) -> NDArray:
149 """
150 Given an array of quad faces return them as triangle faces,
151 also handles pure triangles and mixed triangles and quads.
153 Parameters
154 -----------
155 quads: (n, 4) int
156 Vertex indices of quad faces.
157 dtype
158 Data type requested for the return
159 use_fan
160 Triangulate holes larger than quads with fans,
161 which may be wrong if the holes are non-convex
163 Returns
164 -----------
165 faces : (m, 3) int
166 Vertex indices of triangular faces.
167 """
169 if len(quads) == 0:
170 return np.zeros(0, dtype=dtype)
172 try:
173 # this will fail in newer versions of numpy
174 # if there are mixed quads and tris
175 quads = np.array(quads, dtype=dtype)
177 if len(quads.shape) == 2 and quads.shape[1] == 3:
178 # if they are just triangles return immediately
179 return quads.astype(dtype)
181 if len(quads.shape) == 2 and quads.shape[1] == 4:
182 # if they are just quads stack and return
183 return np.vstack((quads[:, [0, 1, 2]], quads[:, [2, 3, 0]])).astype(dtype)
184 except ValueError:
185 # new numpy raises an error for sequences
186 pass
188 # we made it here so we have mixed tris/quads/polygons
189 # do one pass to get the lengths
190 lengths = np.array([len(i) for i in quads], dtype=np.int64)
192 # get triangles and quads as clean constant-row numpy arrays
193 tri = np.array([quads[i] for i in np.nonzero(lengths == 3)[0]], dtype=np.int64)
194 quad = np.array([quads[i] for i in np.nonzero(lengths == 4)[0]], dtype=np.int64)
196 if use_fan:
197 # get arbitrary polygons as a ragged sequence
198 poly = [quads[i] for i in np.nonzero(lengths > 4)[0]]
199 else:
200 # skip any hole larger than a quad
201 poly = []
203 if len(quad) == 0 and len(poly) == 0:
204 # only triangles, return triangles
205 return tri.astype(dtype)
206 if len(poly) > 0:
207 # use numpy slicing to triangulate ragged-sequence polygons
208 poly = util.triangle_fans_to_faces(poly)
209 if len(quad) > 0:
210 # trivially tessellate quads
211 quad = np.vstack((quad[:, [0, 1, 2]], quad[:, [2, 3, 0]]))
212 # stack triangles from all three cases
213 return util.vstack_empty([tri, quad, poly]).astype(dtype)
216def vertex_face_indices(vertex_count, faces, faces_sparse):
217 """
218 Find vertex face indices from the faces array of vertices
220 Parameters
221 -----------
222 vertex_count : int
223 The number of vertices faces refer to
224 faces : (n, 3) int
225 List of vertex indices
226 faces_sparse : scipy.sparse.COO
227 Sparse matrix
229 Returns
230 -----------
231 vertex_faces : (vertex_count, ) int
232 Face indices for every vertex
233 Array padded with -1 in each row for all vertices with fewer
234 face indices than the max number of face indices.
235 """
236 # Create 2D array with row for each vertex and
237 # length of max number of faces for a vertex
238 try:
239 counts = np.bincount(faces.flatten(), minlength=vertex_count)
240 except TypeError:
241 # casting failed on 32 bit Windows
242 log.warning("casting failed, falling back!")
243 # fall back to np.unique (usually ~35x slower than bincount)
244 counts = np.unique(faces.flatten(), return_counts=True)[1]
245 assert len(counts) == vertex_count
246 assert faces.max() < vertex_count
248 # start cumulative sum at zero and clip off the last value
249 starts = np.append(0, np.cumsum(counts)[:-1])
250 # pack incrementing array into final shape
251 pack = np.arange(counts.max()) + starts[:, None]
252 # pad each row with -1 to pad to the max length
253 padded = -(pack >= (starts + counts)[:, None]).astype(np.int64)
255 try:
256 # do most of the work with a sparse dot product
257 identity = scipy.sparse.identity(len(faces), dtype=int)
258 sorted_faces = faces_sparse.dot(identity).nonzero()[1]
259 # this will fail if any face was degenerate
260 # TODO
261 # figure out how to filter out degenerate faces from sparse
262 # result if sorted_faces.size != faces.size
263 padded[padded == 0] = sorted_faces
264 except BaseException:
265 # fall back to a slow loop
266 log.warning(
267 "vertex_faces falling back to slow loop! "
268 + "mesh probably has degenerate faces",
269 exc_info=True,
270 )
271 sort = np.zeros(faces.size, dtype=np.int64)
272 flat = faces.flatten()
273 for v in range(vertex_count):
274 # assign the data in order
275 sort[starts[v] : starts[v] + counts[v]] = (np.where(flat == v)[0] // 3)[::-1]
276 padded[padded == 0] = sort
277 return padded
280def mean_vertex_normals(vertex_count, faces, face_normals, sparse=None, **kwargs):
281 """
282 Find vertex normals from the mean of the faces that contain
283 that vertex.
285 Parameters
286 -----------
287 vertex_count : int
288 The number of vertices faces refer to
289 faces : (n, 3) int
290 List of vertex indices
291 face_normals : (n, 3) float
292 Normal vector for each face
294 Returns
295 -----------
296 vertex_normals : (vertex_count, 3) float
297 Normals for every vertex
298 Vertices unreferenced by faces will be zero.
299 """
301 def summed_sparse():
302 # use a sparse matrix of which face contains each vertex to
303 # figure out the summed normal at each vertex
304 # allow cached sparse matrix to be passed
305 if sparse is None:
306 matrix = index_sparse(vertex_count, faces)
307 else:
308 matrix = sparse
309 summed = matrix.dot(face_normals)
310 return summed
312 def summed_loop():
313 # loop through every face, in tests was ~50x slower than
314 # doing this with a sparse matrix
315 summed = np.zeros((vertex_count, 3))
316 for face, normal in zip(faces, face_normals):
317 summed[face] += normal
318 return summed
320 try:
321 summed = summed_sparse()
322 except BaseException:
323 log.warning("unable to use sparse matrix, falling back!", exc_info=True)
324 summed = summed_loop()
326 # invalid normals will be returned as zero
327 vertex_normals = util.unitize(summed)
329 return vertex_normals
332def weighted_vertex_normals(
333 vertex_count, faces, face_normals, face_angles, use_loop=False
334):
335 """
336 Compute vertex normals from the faces that contain that vertex.
337 The contribution of a face's normal to a vertex normal is the
338 ratio of the corner-angle in which the vertex is, with respect
339 to the sum of all corner-angles surrounding the vertex.
341 Grit Thuerrner & Charles A. Wuethrich (1998)
342 Computing Vertex Normals from Polygonal Facets,
343 Journal of Graphics Tools, 3:1, 43-46
345 Parameters
346 -----------
347 vertex_count : int
348 The number of vertices faces refer to
349 faces : (n, 3) int
350 List of vertex indices
351 face_normals : (n, 3) float
352 Normal vector for each face
353 face_angles : (n, 3) float
354 Angles at each vertex in the face
356 Returns
357 -----------
358 vertex_normals : (vertex_count, 3) float
359 Normals for every vertex
360 Vertices unreferenced by faces will be zero.
361 """
363 def summed_sparse():
364 # use a sparse matrix of which face contains each vertex to
365 # figure out the summed normal at each vertex
366 # allow cached sparse matrix to be passed
367 # fill the matrix with vertex-corner angles as weights
368 matrix = index_sparse(vertex_count, faces, data=face_angles.ravel())
369 return matrix.dot(face_normals)
371 def summed_loop():
372 summed = np.zeros((vertex_count, 3), np.float64)
373 for vertex_idx in np.arange(vertex_count):
374 # loop over all vertices
375 # compute normal contributions from surrounding faces
376 # obviously slower than with the sparse matrix
377 face_idxs, inface_idxs = np.where(faces == vertex_idx)
378 surrounding_angles = face_angles[face_idxs, inface_idxs]
379 summed[vertex_idx] = np.dot(
380 surrounding_angles / surrounding_angles.sum(), face_normals[face_idxs]
381 )
383 return summed
385 # normals should be unit vectors
386 face_ok = (face_normals**2).sum(axis=1) > 0.5
387 # don't consider faces with invalid normals
388 faces = faces[face_ok]
389 face_normals = face_normals[face_ok]
390 face_angles = face_angles[face_ok]
392 if not use_loop:
393 try:
394 return util.unitize(summed_sparse())
395 except BaseException:
396 log.warning("unable to use sparse matrix, falling back!", exc_info=True)
397 # we either crashed or were asked to loop
398 return util.unitize(summed_loop())
401def index_sparse(columns, indices, data=None, dtype=None):
402 """
403 Return a sparse matrix for which vertices are contained in which faces.
404 A data vector can be passed which is then used instead of booleans
406 Parameters
407 ------------
408 columns : int
409 Number of columns, usually number of vertices
410 indices : (m, d) int
411 Usually mesh.faces
413 Returns
414 ---------
415 sparse: scipy.sparse.coo_matrix of shape (columns, len(faces))
416 dtype is boolean
418 Examples
419 ----------
420 In [1]: sparse = faces_sparse(len(mesh.vertices), mesh.faces)
422 In [2]: sparse.shape
423 Out[2]: (12, 20)
425 In [3]: mesh.faces.shape
426 Out[3]: (20, 3)
428 In [4]: mesh.vertices.shape
429 Out[4]: (12, 3)
431 In [5]: dense = sparse.toarray().astype(int)
433 In [6]: dense
434 Out[6]:
435 array([[1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
436 [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
437 [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
438 [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0],
439 [0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
440 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0],
441 [0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
442 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1],
443 [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0],
444 [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0],
445 [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1],
446 [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1]])
448 In [7]: dense.sum(axis=0)
449 Out[7]: array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3])
450 """
451 indices = np.asanyarray(indices)
452 columns = int(columns)
454 # flattened list
455 row = indices.reshape(-1)
456 col = np.tile(
457 np.arange(len(indices)).reshape((-1, 1)), (1, indices.shape[1])
458 ).reshape(-1)
460 shape = (columns, len(indices))
461 if data is None:
462 data = np.ones(len(col), dtype=bool)
463 elif len(data) != len(col):
464 raise ValueError("data incorrect length")
466 if dtype is not None:
467 data = data.astype(dtype)
469 # assemble into sparse matrix
470 return scipy.sparse.coo_matrix((data, (row, col)), shape=shape, dtype=data.dtype)