Add loop subdivision algorithm

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')

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):
         """

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