Coverage for trimesh/remesh.py: 99%
157 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2remesh.py
3-------------
5Deal with re- triangulation of existing meshes.
6"""
8from itertools import zip_longest
10import numpy as np
12from . import graph, grouping
13from .constants import tol
14from .geometry import faces_to_edges
17def subdivide(
18 vertices, faces, face_index=None, vertex_attributes=None, return_index=False
19):
20 """
21 Subdivide a mesh into smaller triangles.
23 Note that if `face_index` is passed, only those
24 faces will be subdivided and their neighbors won't
25 be modified making the mesh no longer "watertight."
27 Parameters
28 ------------
29 vertices : (n, 3) float
30 Vertices in space
31 faces : (m, 3) int
32 Indexes of vertices which make up triangular faces
33 face_index : faces to subdivide.
34 if None: all faces of mesh will be subdivided
35 if (n,) int array of indices: only specified faces
36 vertex_attributes : dict
37 Contains (n, d) attribute data
38 return_index : bool
39 If True, return index of original face for new faces
41 Returns
42 ----------
43 new_vertices : (q, 3) float
44 Vertices in space
45 new_faces : (p, 3) int
46 Remeshed faces
47 index_dict : dict
48 Only returned if `return_index`, {index of
49 original face : index of new faces}.
50 """
51 if face_index is None:
52 face_mask = np.ones(len(faces), dtype=bool)
53 else:
54 face_mask = np.zeros(len(faces), dtype=bool)
55 face_mask[face_index] = True
57 # the (c, 3) int array of vertex indices
58 faces_subset = faces[face_mask]
60 # find the unique edges of our faces subset
61 edges = np.sort(faces_to_edges(faces_subset), axis=1)
62 unique, inverse = grouping.unique_rows(edges)
63 # then only produce one midpoint per unique edge
64 mid = vertices[edges[unique]].mean(axis=1)
65 mid_idx = inverse.reshape((-1, 3)) + len(vertices)
67 # the new faces_subset with correct winding
68 f = np.column_stack(
69 [
70 faces_subset[:, 0],
71 mid_idx[:, 0],
72 mid_idx[:, 2],
73 mid_idx[:, 0],
74 faces_subset[:, 1],
75 mid_idx[:, 1],
76 mid_idx[:, 2],
77 mid_idx[:, 1],
78 faces_subset[:, 2],
79 mid_idx[:, 0],
80 mid_idx[:, 1],
81 mid_idx[:, 2],
82 ]
83 ).reshape((-1, 3))
85 # add the 3 new faces_subset per old face all on the end
86 # by putting all the new faces after all the old faces
87 # it makes it easier to understand the indexes
88 new_faces = np.vstack((faces[~face_mask], f))
89 # stack the new midpoint vertices on the end
90 new_vertices = np.vstack((vertices, mid))
92 if vertex_attributes is not None:
93 new_attributes = {}
94 for key, values in vertex_attributes.items():
95 if len(values) != len(vertices):
96 continue
97 attr_mid = values[edges[unique]].mean(axis=1)
98 new_attributes[key] = np.vstack((values, attr_mid))
99 return new_vertices, new_faces, new_attributes
101 if return_index:
102 # turn the mask back into integer indexes
103 nonzero = np.nonzero(face_mask)[0]
104 # new faces start past the original faces
105 # but we've removed all the faces in face_mask
106 start = len(faces) - len(nonzero)
107 # indexes are just offset from start
108 stack = np.arange(start, start + len(f) * 4).reshape((-1, 4))
109 # reformat into a slightly silly dict for some reason
110 index_dict = dict(zip(nonzero, stack))
112 return new_vertices, new_faces, index_dict
114 return new_vertices, new_faces
117def subdivide_to_size(vertices, faces, max_edge, max_iter=10, return_index=False):
118 """
119 Subdivide a mesh until every edge is shorter than a
120 specified length.
122 Every edge longer than `max_edge` is bisected at a single shared
123 midpoint, so the two faces on either side stay in sync and a
124 watertight input stays watertight — no T-junctions (cracks) are
125 introduced. Faces already small enough are left untouched.
127 Parameters
128 ------------
129 vertices : (n, 3) float
130 Vertices in space
131 faces : (m, 3) int
132 Indices of vertices which make up triangles
133 max_edge : float
134 Maximum length of any edge in the result
135 max_iter : int
136 The maximum number of times to run subdivision
137 return_index : bool
138 If True, return index of original face for new faces
140 Returns
141 ------------
142 vertices : (j, 3) float
143 Vertices in space
144 faces : (q, 3) int
145 Indices of vertices
146 index : (q,) int
147 Only returned if `return_index`, index of
148 original face for each new face.
149 """
150 # copy inputs and make sure dtype is correct
151 vertices = np.array(vertices, dtype=np.float64, copy=True)
152 faces = np.array(faces, dtype=np.int64, copy=True)
153 # map each current face back to its original face index
154 index = np.arange(len(faces))
156 for i in range(max_iter + 1):
157 # the unique edges of the current mesh and the length of each
158 edges = np.sort(faces_to_edges(faces), axis=1)
159 unique, inverse = grouping.unique_rows(edges)
160 edges_unique = edges[unique]
161 # length uses xyz only — callers may hstack extra columns (e.g. uv)
162 lengths = (
163 (vertices[edges_unique[:, 0], :3] - vertices[edges_unique[:, 1], :3]) ** 2
164 ).sum(axis=1) ** 0.5
165 long_edge = lengths > max_edge
167 # every edge is short enough so we're done
168 if not long_edge.any():
169 break
170 # check max_iter before refining again
171 if i >= max_iter:
172 raise ValueError("max_iter exceeded!")
174 # one shared midpoint vertex per long edge, -1 where left intact
175 midpoint = np.full(len(unique), -1, dtype=np.int64)
176 midpoint[long_edge] = np.arange(int(long_edge.sum())) + len(vertices)
177 vertices = np.vstack((vertices, vertices[edges_unique[long_edge]].mean(axis=1)))
179 # midpoint id for each face's 3 edges, columns (v0v1, v1v2, v2v0)
180 face_mid = midpoint[inverse].reshape((-1, 3))
181 split = face_mid >= 0
182 # number of edges marked for bisection on each face: 0, 1, 2, or 3
183 count = split.sum(axis=1)
185 # 0 marked edges — face passes through unchanged
186 keep = count == 0
187 faces_keep = faces[keep]
188 index_keep = index[keep]
190 # 1 marked edge — rotate so the split edge is (a, b) so a single
191 # template handles all three cases: fan the midpoint to the
192 # opposite vertex as [a, p, c], [p, b, c]
193 one = count == 1
194 j = np.argmax(split[one], axis=1)
195 r = np.arange(len(j))
196 f1, mid1 = faces[one], face_mid[one]
197 a, b, c = f1[r, j], f1[r, (j + 1) % 3], f1[r, (j + 2) % 3]
198 p = mid1[r, j]
199 faces_one = np.column_stack((a, p, c, p, b, c)).reshape((-1, 3))
200 index_one = np.repeat(index[one], 2)
202 # 2 marked edges — rotate so the unsplit edge is (c, a), again so
203 # one template covers all three cases: a corner triangle [p, b, q]
204 # plus the quad (a, p, q, c) cut along its shorter diagonal
205 two = count == 2
206 j = (np.argmin(split[two], axis=1) + 1) % 3
207 r = np.arange(len(j))
208 f2, mid2 = faces[two], face_mid[two]
209 a, b, c = f2[r, j], f2[r, (j + 1) % 3], f2[r, (j + 2) % 3]
210 p, q = mid2[r, j], mid2[r, (j + 1) % 3]
211 # the shorter of the two quad diagonals: a-q versus p-c (xyz only)
212 use_aq = (
213 ((vertices[a, :3] - vertices[q, :3]) ** 2).sum(axis=1)
214 <= ((vertices[p, :3] - vertices[c, :3]) ** 2).sum(axis=1)
215 )[:, None]
216 corner = np.column_stack((p, b, q))
217 quad_a = np.where(use_aq, np.column_stack((a, p, q)), np.column_stack((a, p, c)))
218 quad_b = np.where(use_aq, np.column_stack((a, q, c)), np.column_stack((p, q, c)))
219 faces_two = np.vstack((corner, quad_a, quad_b))
220 index_two = np.tile(index[two], 3)
222 # 3 marked edges — the regular 1 -> 4 split, same winding as subdivide
223 three = count == 3
224 f3, mid3 = faces[three], face_mid[three]
225 v0, v1, v2 = f3.T
226 m0, m1, m2 = mid3.T
227 faces_three = np.column_stack(
228 (v0, m0, m2, m0, v1, m1, m2, m1, v2, m0, m1, m2)
229 ).reshape((-1, 3))
230 index_three = np.repeat(index[three], 4)
232 faces = np.vstack((faces_keep, faces_one, faces_two, faces_three)).astype(
233 np.int64
234 )
235 index = np.concatenate((index_keep, index_one, index_two, index_three))
237 if return_index:
238 assert len(index) == len(faces)
239 return vertices, faces, index
241 return vertices, faces
244def subdivide_loop(vertices, faces, iterations=None):
245 """
246 Subdivide a mesh by dividing each triangle into four triangles
247 and approximating their smoothed surface (loop subdivision).
248 This function is an array-based implementation of loop subdivision,
249 which avoids slow for loop and enables faster calculation.
251 Overall process:
252 1. Calculate odd vertices.
253 Assign a new odd vertex on each edge and
254 calculate the value for the boundary case and the interior case.
255 The value is calculated as follows.
256 v2
257 / f0 \\ 0
258 v0--e--v1 / \\
259 \\f1 / v0--e--v1
260 v3
261 - interior case : 3:1 ratio of mean(v0,v1) and mean(v2,v3)
262 - boundary case : mean(v0,v1)
263 2. Calculate even vertices.
264 The new even vertices are calculated with the existing
265 vertices and their adjacent vertices.
266 1---2
267 / \\/ \\ 0---1
268 0---v---3 / \\/ \\
269 \\ /\\/ b0---v---b1
270 k...4
271 - interior case : (1-kB):B ratio of v and k adjacencies
272 - boundary case : 3:1 ratio of v and mean(b0,b1)
273 3. Compose new faces with new vertices.
275 Parameters
276 ------------
277 vertices : (n, 3) float
278 Vertices in space
279 faces : (m, 3) int
280 Indices of vertices which make up triangles
282 Returns
283 ------------
284 vertices : (j, 3) float
285 Vertices in space
286 faces : (q, 3) int
287 Indices of vertices
288 iterations : int
289 Number of iterations to run subdivision
290 """
291 if iterations is None:
292 iterations = 1
294 def _subdivide(vertices, faces):
295 # find the unique edges of our faces
296 edges, edges_face = faces_to_edges(faces, return_index=True)
297 edges.sort(axis=1)
298 unique, inverse = grouping.unique_rows(edges)
300 # set interior edges if there are two edges and boundary if there is
301 # one.
302 edge_inter = np.sort(grouping.group_rows(edges, require_count=2), axis=1)
303 edge_bound = grouping.group_rows(edges, require_count=1)
304 # make sure that one edge is shared by only one or two faces.
305 if not len(edge_inter) * 2 + len(edge_bound) == len(edges):
306 # we have multiple bodies it's a party!
307 # edges shared by 2 faces are "connected"
308 # so this connected components operation is
309 # essentially identical to `face_adjacency`
310 faces_group = graph.connected_components(edges_face[edge_inter])
312 if len(faces_group) == 1:
313 raise ValueError("Some edges are shared by more than 2 faces")
315 # collect a subdivided copy of each body
316 seq_verts = []
317 seq_faces = []
318 # keep track of vertex count as we go so
319 # we can do a single vstack at the end
320 count = 0
321 # loop through original face indexes
322 for f in faces_group:
323 # a lot of the complexity in this operation
324 # is computing vertex neighbors so we only
325 # want to pass forward the referenced vertices
326 # for this particular group of connected faces
327 unique, inverse = grouping.unique_bincount(
328 faces[f].reshape(-1), return_inverse=True
329 )
331 # subdivide this subset of faces
332 cur_verts, cur_faces = _subdivide(
333 vertices=vertices[unique], faces=inverse.reshape((-1, 3))
334 )
336 # increment the face references to match
337 # the vertices when we stack them later
338 cur_faces += count
339 # increment the total vertex count
340 count += len(cur_verts)
341 # append to the sequence
342 seq_verts.append(cur_verts)
343 seq_faces.append(cur_faces)
345 # return results as clean (n, 3) arrays
346 return np.vstack(seq_verts), np.vstack(seq_faces)
348 # set interior, boundary mask for unique edges
349 edge_bound_mask = np.zeros(len(edges), dtype=bool)
350 edge_bound_mask[edge_bound] = True
351 edge_bound_mask = edge_bound_mask[unique]
352 edge_inter_mask = ~edge_bound_mask
354 # find the opposite face for each edge
355 edge_pair = np.zeros(len(edges)).astype(int)
356 edge_pair[edge_inter[:, 0]] = edge_inter[:, 1]
357 edge_pair[edge_inter[:, 1]] = edge_inter[:, 0]
358 opposite_face1 = edges_face[unique]
359 opposite_face2 = edges_face[edge_pair[unique]]
361 # set odd vertices to the middle of each edge (default as boundary
362 # case).
363 odd = vertices[edges[unique]].mean(axis=1)
364 # modify the odd vertices for the interior case
365 e = edges[unique[edge_inter_mask]]
366 e_v0 = vertices[e][:, 0]
367 e_v1 = vertices[e][:, 1]
368 e_f0 = faces[opposite_face1[edge_inter_mask]]
369 e_f1 = faces[opposite_face2[edge_inter_mask]]
370 e_v2_idx = e_f0[~(e_f0[:, :, None] == e[:, None, :]).any(-1)]
371 e_v3_idx = e_f1[~(e_f1[:, :, None] == e[:, None, :]).any(-1)]
372 e_v2 = vertices[e_v2_idx]
373 e_v3 = vertices[e_v3_idx]
375 # simplified from:
376 # # 3 / 8 * (e_v0 + e_v1) + 1 / 8 * (e_v2 + e_v3)
377 odd[edge_inter_mask] = 0.375 * e_v0 + 0.375 * e_v1 + e_v2 / 8.0 + e_v3 / 8.0
379 # find vertex neighbors of each vertex
380 neighbors = graph.neighbors(edges=edges[unique], max_index=len(vertices))
381 # convert list type of array into a fixed-shaped numpy array (set -1 to
382 # empties)
383 neighbors = np.array(list(zip_longest(*neighbors, fillvalue=-1))).T
384 # if the neighbor has -1 index, its point is (0, 0, 0), so that
385 # it is not included in the summation of neighbors when calculating the
386 # even
387 vertices_ = np.vstack([vertices, [0.0, 0.0, 0.0]])
388 # number of neighbors
389 k = (neighbors + 1).astype(bool).sum(axis=1)
391 # calculate even vertices for the interior case
392 even = np.zeros_like(vertices)
394 # beta = 1 / k * (5 / 8 - (3 / 8 + 1 / 4 * np.cos(2 * np.pi / k)) ** 2)
395 # simplified with sympy.parse_expr('...').simplify()
396 beta = (40.0 - (2.0 * np.cos(2 * np.pi / k) + 3) ** 2) / (64 * k)
397 even = (
398 beta[:, None] * vertices_[neighbors].sum(1)
399 + (1 - k[:, None] * beta[:, None]) * vertices
400 )
402 # calculate even vertices for the boundary case
403 if edge_bound_mask.any():
404 # boundary vertices from boundary edges
405 vrt_bound_mask = np.zeros(len(vertices), dtype=bool)
406 vrt_bound_mask[np.unique(edges[unique][~edge_inter_mask])] = True
407 # one boundary vertex has two neighbor boundary vertices (set
408 # others as -1)
409 boundary_neighbors = neighbors[vrt_bound_mask]
410 boundary_neighbors[~vrt_bound_mask[neighbors[vrt_bound_mask]]] = -1
412 even[vrt_bound_mask] = (
413 vertices_[boundary_neighbors].sum(axis=1) / 8.0
414 + (3.0 / 4.0) * vertices[vrt_bound_mask]
415 )
417 # the new faces with odd vertices
418 odd_idx = inverse.reshape((-1, 3)) + len(vertices)
419 new_faces = np.column_stack(
420 [
421 faces[:, 0],
422 odd_idx[:, 0],
423 odd_idx[:, 2],
424 odd_idx[:, 0],
425 faces[:, 1],
426 odd_idx[:, 1],
427 odd_idx[:, 2],
428 odd_idx[:, 1],
429 faces[:, 2],
430 odd_idx[:, 0],
431 odd_idx[:, 1],
432 odd_idx[:, 2],
433 ]
434 ).reshape((-1, 3))
436 # stack the new even vertices and odd vertices
437 new_vertices = np.vstack((even, odd))
439 return new_vertices, new_faces
441 for _ in range(iterations):
442 vertices, faces = _subdivide(vertices, faces)
444 if tol.strict or True:
445 assert np.isfinite(vertices).all()
446 assert np.isfinite(faces).all()
447 # should raise if faces are malformed
448 assert np.isfinite(vertices[faces]).all()
450 # none of the faces returned should be degenerate
451 # i.e. every face should have 3 unique vertices
452 assert (faces[:, 1:] != faces[:, :1]).all()
454 return vertices, faces