diff --git a/tests/test_remesh.py b/tests/test_remesh.py index 6411e338a..7e0848d2d 100644 --- a/tests/test_remesh.py +++ b/tests/test_remesh.py @@ -106,6 +106,56 @@ def test_sub(self): # volume should be the same assert g.np.isclose(m.volume, s.volume) + def test_loop(self): + meshes = [ + g.get_mesh('soup.stl'), # a soup of random triangles + g.get_mesh('featuretype.STL')] # a mesh with a single body + + for m in meshes: + sub = m.loop(iterations=1) + # number of faces should increase + assert len(sub.faces) > len(m.faces) + # subdivided faces are smaller + assert sub.area_faces.mean() < m.area_faces.mean() + + def test_loop_multibody(self): + mesh = g.get_mesh('cycloidal.ply') # a mesh with multiple bodies + sub = mesh.loop(iterations=1, multibody=True) + # number of faces should increase + assert len(sub.faces) > len(mesh.faces) + # subdivided faces are smaller + assert sub.area_faces.mean() < mesh.area_faces.mean() + + def test_loop_correct(self): + box = g.trimesh.creation.box() + big_sphere = g.trimesh.creation.icosphere(radius=0.5) + small_sphere = g.trimesh.creation.icosphere(radius=0.4) + sub = box.loop(iterations=2) + # smaller than 0.5 sphere + assert big_sphere.contains(sub.vertices).all() + # bigger than 0.4 sphere + assert (~small_sphere.contains(sub.vertices)).all() + + def test_loop_bound(self): + def _get_boundary_vertices(mesh): + boundary_groups = g.trimesh.grouping.group_rows( + mesh.edges_sorted, require_count=1) + return mesh.vertices[g.np.unique(mesh.edges_sorted[boundary_groups])] + + box = g.trimesh.creation.box() + bottom_mask = g.np.zeros(len(box.faces), dtype=bool) + bottom_faces = [1, 5] + bottom_mask[bottom_faces] = True + # eliminate bottom of the box + box.update_faces(~bottom_mask) + bottom_vrts = _get_boundary_vertices(box) + # subdivide box + sub = box.loop(iterations=2) + sub_bottom_vrts = _get_boundary_vertices(sub) + epsilon = 1e-5 + # y value of bottom boundary vertices should not be changed + assert (bottom_vrts[:, 1].mean() - sub_bottom_vrts[:, 1].mean()) < epsilon + def test_uv(self): # get a mesh with texture m = g.get_mesh('fuze.obj') diff --git a/trimesh/base.py b/trimesh/base.py index 5e63e6bc0..bd101225f 100644 --- a/trimesh/base.py +++ b/trimesh/base.py @@ -1997,6 +1997,48 @@ def subdivide_to_size(self, max_edge, max_iter=10, return_index=False): return result + def loop(self, iterations=1, multibody=False): + """ + Subdivide a mesh by dividing each triangle into four triangles + and approximating their smoothed surface (loop subdivision). + + Parameters + ------------ + iterations : int + Number of iterations to run subdivisio + multibody : bool + If True will try to subdivide for each submesh + """ + if multibody: + splited_meshes = self.split(only_watertight=False) + if len(splited_meshes) > 1: + new_meshes = [] + # perform subdivision for all submesh + for splited_mesh in splited_meshes: + new_vertices, new_faces = remesh.loop( + vertices=splited_mesh.vertices, + faces=splited_mesh.faces, + iterations=iterations) + # create new mesh + new_mesh = Trimesh( + vertices=new_vertices, + faces=new_faces) + new_meshes.append(new_mesh) + # concatenate all meshes into one + result = util.concatenate(new_meshes) + return result + + # perform subdivision for one mesh + new_vertices, new_faces = remesh.loop( + vertices=self.vertices, + faces=self.faces, + iterations=iterations) + # create new mesh + result = Trimesh( + vertices=new_vertices, + faces=new_faces) + return result + @log_time def smoothed(self, **kwargs): """ diff --git a/trimesh/remesh.py b/trimesh/remesh.py index d162eb4ff..60bcb8234 100644 --- a/trimesh/remesh.py +++ b/trimesh/remesh.py @@ -4,11 +4,11 @@ Deal with re- triangulation of existing meshes. """ - import numpy as np from . import util from . import grouping +from . import graph from .geometry import faces_to_edges @@ -213,3 +213,150 @@ def subdivide_to_size(vertices, return final_vertices, final_faces, final_index return final_vertices, final_faces + + +def loop(vertices, + faces, + iterations=1): + """ + Subdivide a mesh by dividing each triangle into four triangles + and approximating their smoothed surface (loop subdivision). + This function is an array-based implementation of loop subdivision, + which avoids slow for loop and enables faster calculation. + + Overall process: + 1. Calculate odd vertices. + Assign a new odd vertex on each edge and + calculate the value for the boundary case and the interior case. + The value is calculated as follows. + v2 + / f0 \\ 0 + v0--e--v1 / \\ + \\f1 / v0--e--v1 + v3 + - interior case : 3:1 ratio of mean(v0,v1) and mean(v2,v3) + - boundary case : mean(v0,v1) + 2. Calculate even vertices. + The new even vertices are calculated with the existing + vertices and their adjacent vertices. + 1---2 + / \\/ \\ 0---1 + 0---v---3 / \\/ \\ + \\ /\\/ b0---v---b1 + k...4 + - interior case : (1-kβ):β ratio of v and k adjacencies + - boundary case : 3:1 ratio of v and mean(b0,b1) + 3. Compose new faces with new vertices. + + Parameters + ------------ + vertices : (n, 3) float + Vertices in space + faces : (m, 3) int + Indices of vertices which make up triangles + + Returns + ------------ + vertices : (j, 3) float + Vertices in space + faces : (q, 3) int + Indices of vertices + iterations : int + Number of iterations to run subdivision + """ + try: + from itertools import zip_longest + except BaseException: + from itertools import izip_longest as zip_longest # python2 + + def _subdivide(vertices, faces): + # find the unique edges of our faces + edges, edges_face = faces_to_edges(faces, return_index=True) + edges = np.sort(edges, axis=1) + unique, inverse = grouping.unique_rows(edges) + + # set interior edges if there are two edges and boundary if there is one. + edge_inter = np.sort(grouping.group_rows(edges, require_count=2), axis=1) + edge_bound = grouping.group_rows(edges, require_count=1) + # make sure that one edge is shared by only one or two faces. + if not len(edge_inter)*2 + len(edge_bound) == len(edges): + raise ValueError('Some edges are shared by more than 2 faces') + + # set interior, boundary mask for unique edges + edge_bound_mask = np.zeros(len(edges), dtype=bool) + edge_bound_mask[edge_bound] = True + edge_bound_mask = edge_bound_mask[unique] + edge_inter_mask = ~edge_bound_mask + + # find the opposite face for each edge + edge_pair = np.zeros(len(edges)).astype(int) + edge_pair[edge_inter[:, 0]] = edge_inter[:, 1] + edge_pair[edge_inter[:, 1]] = edge_inter[:, 0] + opposite_face1 = edges_face[unique] + opposite_face2 = edges_face[edge_pair[unique]] + + # set odd vertices to the middle of each edge (default as boundary case). + odd = vertices[edges[unique]].mean(axis=1) + # modify the odd vertices for the interior case + e = edges[unique[edge_inter_mask]] + e_v0 = vertices[e][:, 0] + e_v1 = vertices[e][:, 1] + e_f0 = faces[opposite_face1[edge_inter_mask]] + e_f1 = faces[opposite_face2[edge_inter_mask]] + e_v2_idx = e_f0[~(e_f0[:, :, None] == e[:, None, :]).any(-1)] + e_v3_idx = e_f1[~(e_f1[:, :, None] == e[:, None, :]).any(-1)] + e_v2 = vertices[e_v2_idx] + e_v3 = vertices[e_v3_idx] + odd[edge_inter_mask] = 3/8 * (e_v0 + e_v1) + 1/8 * (e_v2 + e_v3) + + # find vertex neighbors of each vertex + neighbors = graph.neighbors(edges=edges[unique], max_index=len(vertices)) + # convert list type of array into a fixed-shaped numpy array (set -1 to empties) + neighbors = np.array(list(zip_longest(*neighbors, fillvalue=-1))).T + # if the neighbor has -1 index, its point is (0, 0, 0), so that + # it is not included in the summation of neighbors when calculating the even + vertices_ = np.vstack([vertices, [0, 0, 0]]) + # number of neighbors + k = (neighbors + 1).astype(bool).sum(-1) + + # calculate even vertices for the interior case + even = np.zeros_like(vertices) + beta = 1/k * (5/8 - (3/8 + 1/4 * np.cos(2 * np.pi / k)) ** 2) + even = beta[:, None] * vertices_[neighbors].sum(1) \ + + (1 - k[:, None] * beta[:, None]) * vertices + + # calculate even vertices for the boundary case + if True in edge_bound_mask: + # boundary vertices from boundary edges + vrt_bound_mask = np.zeros(len(vertices), dtype=bool) + vrt_bound_mask[np.unique(edges[unique][~edge_inter_mask])] = True + # one boundary vertex has two neighbor boundary vertices (set others as -1) + boundary_neighbors = neighbors[vrt_bound_mask] + boundary_neighbors[~vrt_bound_mask[neighbors[vrt_bound_mask]]] = -1 + even[vrt_bound_mask] = 1/8 * vertices_[boundary_neighbors].sum(1) \ + + 3/4 * vertices[vrt_bound_mask] + + # the new faces with odd vertices + odd_idx = inverse.reshape((-1, 3)) + len(vertices) + new_faces = np.column_stack([ + faces[:, 0], + odd_idx[:, 0], + odd_idx[:, 2], + odd_idx[:, 0], + faces[:, 1], + odd_idx[:, 1], + odd_idx[:, 2], + odd_idx[:, 1], + faces[:, 2], + odd_idx[:, 0], + odd_idx[:, 1], + odd_idx[:, 2]]).reshape((-1, 3)) + + # stack the new even vertices and odd vertices + new_vertices = np.vstack((even, odd)) + + return new_vertices, new_faces + + for _index in range(iterations): + vertices, faces = _subdivide(vertices, faces) + return vertices, faces