Repository: mikedh/trimesh
Pull request #265 (2018-11-30 22:12:49)

Slice with multiple planes and performance tweaks
Two main things:
1) I found that slicing with multiple planes (not necessarily parallel) gets a speed up when not creating a new Trimesh object every time, so I've expanded the slice_plane method to slice_planes, where you can instead pass an array of origins/normals. I left the old method for compatibility sake, but feel free to do what you want with it (e.g., the new method subsumes the old, so eventually might want to just switch over to the new one)!
2) Some minor performance tweaks. The biggest one was replacing np.unique with np.bincount. Since the faces arrays are all integer-valued, it seems like np.bincount/cumsum perform a lot better than np.unique (e.g., 2-3x faster for small meshes, 10x faster for large meshes). Here are some examples for a smaller mesh (10s of faces) and a larger mesh (~40,000) faces (using kernprof):
![headings](https://user-images.githubusercontent.com/7028710/49317488-8999ab80-f4a9-11e8-9521-4393c7b938fa.png)
![np_unique](https://user-images.githubusercontent.com/7028710/49317492-8d2d3280-f4a9-11e8-839e-2ac8b2943026.png)
![np_unique_big](https://user-images.githubusercontent.com/7028710/49317496-8ef6f600-f4a9-11e8-9368-0a15c807ac64.png)
![unique_replacement](https://user-images.githubusercontent.com/7028710/49317499-90c0b980-f4a9-11e8-8e71-e348dd9f625b.png)
![unique_replacement_big](https://user-images.githubusercontent.com/7028710/49317507-9d451200-f4a9-11e8-9954-2a67e61ebce9.png)

Happy to answer questions and I hope it helps!



Patch:
diff --git a/trimesh/base.py b/trimesh/base.py
index 8bdd6f7f..34eaf924 100644
--- a/trimesh/base.py
+++ b/trimesh/base.py
@@ -1858,6 +1858,39 @@ class Trimesh(Geometry):
             **kwargs)
 
         return new_mesh
+    
+    def slice_planes(self,
+                     plane_origins,
+                     plane_normals,
+                     **kwargs):
+        """
+        Returns another mesh that is the current mesh
+        sliced by the planes defined by origins and normals.
+        Saves some computation over calling the slice_mesh_plane
+        method multiple times by not creating a new object and 
+        just carrying the faces and vertices arrays.
+
+        Parameters
+        ---------
+        plane_normals: (n,3) float
+          Normal vectors of slicing planes
+        plane_origins : (n,3) float
+          Points on the slicing planes
+
+        Returns
+        ---------
+        new_mesh: trimesh.Trimesh or None
+          Subset of current mesh sliced by planes
+        """
+
+        # return a new mesh
+        new_mesh = intersections.slice_mesh_planes(
+            mesh=self,
+            plane_normals=plane_normals,
+            plane_origins=plane_origins,
+            **kwargs)
+
+        return new_mesh
 
     @caching.cache_decorator
     def convex_hull(self):
diff --git a/trimesh/intersections.py b/trimesh/intersections.py
index bb756748..e0589df4 100644
--- a/trimesh/intersections.py
+++ b/trimesh/intersections.py
@@ -385,19 +385,22 @@ def planes_lines(plane_origins,
 
     return on_plane, valid
 
-
-def slice_mesh_plane(mesh,
-                     plane_normal,
-                     plane_origin,
-                     cached_dots=None):
+def slice_faces_plane(vertices, 
+                      faces,
+                      plane_normal,
+                      plane_origin,
+                      cached_dots=None):
     """
-    Slice a mesh with a plane, returning a new mesh that is the
-    portion of the original mesh to the positive normal side of the plane
+    Slice a mesh (given as a set of faces and vertices) with a plane, returning a 
+    new mesh (again as a set of faces and vertices) that is the
+    portion of the original mesh to the positive normal side of the plane.
 
     Parameters
     ---------
-    mesh : Trimesh object
-        Source mesh to slice
+    vertices : (n, 3) float
+        Vertices of source mesh to slice
+    faces : (n, 3) int
+        Faces of source mesh to slice
     plane_normal : (3,) float
         Normal vector of plane to intersect with mesh
     plane_origin:  (3,) float
@@ -407,71 +410,68 @@ def slice_mesh_plane(mesh,
         products pass them here to avoid recomputing
     Returns
     ----------
-    new_mesh : Trimesh object
-        Sliced mesh
+    new_vertices : (n, 3) float
+        Vertices of sliced mesh
+    new_faces : (n, 3) int
+        Faces of sliced mesh
     """
 
-    from . import base
-
-    # If passed invalid input, return None
-    if mesh is None:
-        return None
-
-    # check input plane
-    plane_normal = np.asanyarray(plane_normal,
-                                 dtype=np.float64)
-    plane_origin = np.asanyarray(plane_origin,
-                                 dtype=np.float64)
-
-    if plane_origin.shape != (3,) or plane_normal.shape != (3,):
-        raise ValueError('Plane origin and normal must be (3,)!')
+    if len(vertices) == 0:
+        return vertices, faces
 
     if cached_dots is not None:
         dots = cached_dots
     else:
         # dot product of each vertex with the plane normal indexed by face
         # so for each face the dot product of each vertex is a row
-        # shape is the same as mesh.faces (n,3)
-        dots = np.dot(plane_normal,
-                      (mesh.vertices - plane_origin).T)[mesh.faces]
+        # shape is the same as faces (n,3)
+        dots = np.einsum('i,ij->j', plane_normal, 
+                                    (vertices - plane_origin).T)[faces]
 
     # Find vertex orientations w.r.t. faces for all triangles:
     #  -1 -> vertex "inside" plane (positive normal direction)
     #   0 -> vertex on plane
     #   1 -> vertex "outside" plane (negative normal direction)
-    signs = np.zeros(mesh.faces.shape, dtype=np.int8)
+    signs = np.zeros(faces.shape, dtype=np.int8)
     signs[dots < -tol.merge] = 1
     signs[dots > tol.merge] = -1
+    signs[np.logical_and(dots >= -tol.merge, dots <= tol.merge)] = 0
 
     # Find all triangles that intersect this plane
     # onedge <- indices of all triangles intersecting the plane
     # inside <- indices of all triangles "inside" the plane (positive normal)
-    signs_sum = signs.sum(axis=1)
-    signs_asum = np.abs(signs).sum(axis=1)
+    signs_sum = signs.sum(axis=1, dtype=np.int8)
+    signs_asum = np.abs(signs).sum(axis=1, dtype=np.int8)
 
+    # Cases:
+    # (0,0,0),  (-1,0,0),  (-1,-1,0), (-1,-1,-1) <- inside
+    # (1,0,0),  (1,1,0),   (1,1,1)               <- outside
+    # (1,0,-1), (1,-1,-1), (1,1,-1)              <- onedge
     onedge = np.logical_and(signs_asum >= 2,
                             np.abs(signs_sum) <= 1)
-    inside = np.logical_or(signs_sum == -3,
-                           np.logical_and(signs_sum == -1, signs_asum == 1))
+    inside = (signs_sum == -signs_asum)
 
     # Automatically include all faces that are "inside"
-    new_faces = mesh.faces[inside]
+    new_faces = faces[inside]
 
     # Separate faces on the edge into two cases: those which will become
     # quads (two vertices inside plane) and those which will become triangles
     # (one vertex inside plane)
-    cut_triangles = mesh.triangles[onedge]
-    cut_faces_quad = mesh.faces[np.logical_and(onedge, signs_sum < 0)]
-    cut_faces_tri = mesh.faces[np.logical_and(onedge, signs_sum >= 0)]
+    triangles = vertices[faces]
+    cut_triangles = triangles[onedge]
+    cut_faces_quad = faces[np.logical_and(onedge, signs_sum < 0)]
+    cut_faces_tri = faces[np.logical_and(onedge, signs_sum >= 0)]
     cut_signs_quad = signs[np.logical_and(onedge, signs_sum < 0)]
     cut_signs_tri = signs[np.logical_and(onedge, signs_sum >= 0)]
 
     # If no faces to cut, the surface is not in contact with this plane.
     # Thus, return a mesh with only the inside faces
     if len(cut_faces_quad) + len(cut_faces_tri) == 0:
-        new_mesh = base.Trimesh(vertices=mesh.vertices, faces=new_faces)
-        new_mesh.remove_unreferenced_vertices()
-        return new_mesh
+        vert_counts = np.bincount(new_faces.flatten(), minlength=len(vertices))
+        unique_verts = vert_counts > 0
+        unique_inds = np.where(unique_verts)[0]
+        unique_faces = (np.cumsum(unique_verts)-1)[new_faces]
+        return vertices[unique_inds], unique_faces
 
     # Extract the intersections of each triangle's edges with the plane
     o = cut_triangles                               # origins
@@ -484,7 +484,7 @@ def slice_mesh_plane(mesh,
     int_points = np.einsum('ij,ijk->ijk', dist, d) + o
 
     # Initialize the array of new vertices with the current vertices
-    new_vertices = mesh.vertices
+    new_vertices = vertices
 
     # Handle the case where a new quad is formed by the intersection
     # First, extract the intersection points belonging to a new quad
@@ -550,10 +550,103 @@ def slice_mesh_plane(mesh,
         new_faces = np.append(new_faces, new_tri_faces, axis=0)
 
     # new mesh without cut off vertices
-    unique, inverse = np.unique(new_faces,
-                                return_inverse=True)
-    new_mesh = base.Trimesh(vertices=new_vertices[unique],
-                            faces=inverse.reshape((-1, 3)),
+    vert_counts = np.bincount(new_faces.flatten(), minlength=len(new_vertices))
+    unique_verts = vert_counts > 0
+    unique_inds = np.where(unique_verts)[0]
+    unique_faces = (np.cumsum(unique_verts)-1)[new_faces]
+    return new_vertices[unique_inds], unique_faces
+    
+def slice_mesh_plane(mesh,
+                     plane_normal,
+                     plane_origin,
+                     cached_dots=None):
+    """
+    Slice a mesh with a plane, returning a new mesh that is the
+    portion of the original mesh to the positive normal side of the plane
+
+    Parameters
+    ---------
+    mesh : Trimesh object
+        Source mesh to slice
+    plane_normal : (3,) float
+        Normal vector of plane to intersect with mesh
+    plane_origin:  (3,) float
+        Point on plane to intersect with mesh
+    cached_dots : (n, 3) float
+        If an external function has stored dot
+        products pass them here to avoid recomputing
+    Returns
+    ----------
+    new_mesh : Trimesh object
+        Sliced mesh
+    """
+
+    from . import base
+
+    # If passed invalid input, return None
+    if mesh is None:
+        return None
+
+    # check input plane
+    plane_normal = np.asanyarray(plane_normal,
+                                 dtype=np.float64)
+    plane_origin = np.asanyarray(plane_origin,
+                                 dtype=np.float64)
+
+    if plane_origin.shape != (3,) or plane_normal.shape != (3,):
+        raise ValueError('Plane origin and normal must be (3,)!')
+    
+    new_vertices, new_faces = slice_faces_plane(mesh.vertices, mesh.faces, plane_normal, plane_origin, cached_dots=cached_dots)
+    new_mesh = base.Trimesh(vertices=new_vertices,
+                            faces=new_faces,
                             process=False)
+    return new_mesh
 
+def slice_mesh_planes(mesh,
+                      plane_normals,
+                      plane_origins):
+    """
+    Slice a mesh with a plane, returning a new mesh that is the
+    portion of the original mesh to the positive normal side of the plane
+
+    Parameters
+    ---------
+    mesh : Trimesh object
+        Source mesh to slice
+    plane_normal : (3,) float
+        Normal vector of plane to intersect with mesh
+    plane_origin:  (3,) float
+        Point on plane to intersect with mesh
+    cached_dots : (n, 3) float
+        If an external function has stored dot
+        products pass them here to avoid recomputing
+    Returns
+    ----------
+    new_mesh : Trimesh object
+        Sliced mesh
+    """
+
+    from . import base
+
+    # If passed invalid input, return None
+    if mesh is None:
+        return None
+
+    # check input plane
+    plane_normals = np.asanyarray(plane_normals,
+                                dtype=np.float64)
+    plane_origins = np.asanyarray(plane_origins,
+                                dtype=np.float64)
+
+    if plane_origins.shape[1] != 3 or plane_normals.shape[1] != 3:
+        raise ValueError('Plane origins and normals must be of dimension (n,3)!')
+
+    new_vertices = mesh.vertices
+    new_faces = mesh.faces
+    for o,n in zip(plane_origins, plane_normals):
+        new_vertices, new_faces = slice_faces_plane(new_vertices, new_faces, n, o)
+    
+    new_mesh = base.Trimesh(vertices=new_vertices,
+                            faces=new_faces,
+                            process=False)
     return new_mesh


Test patch:
diff --git a/tests/test_section.py b/tests/test_section.py
index 2c8c6391..7f353ff1 100644
--- a/tests/test_section.py
+++ b/tests/test_section.py
@@ -176,6 +176,32 @@ class SliceTest(g.unittest.TestCase):
         # should be lots of stuff at the plane and nothing behind
         assert g.np.isclose(dot.min(), 0.0)
 
+        # Cut part of top off of box with multiple planes at once
+        # and make sure bounds and number of faces is correct
+        plane_origins = [mesh.bounds[1] - 0.05, mesh.bounds[1] - 0.05]
+        plane_normals = [g.np.array([0, 0, 1]), g.np.array([0, 1, 0])]
+
+        sliced = mesh.slice_planes(plane_origins=plane_origins,
+                                   plane_normals=plane_normals)
+
+        assert g.np.isclose(
+            sliced.bounds[0], mesh.bounds[0] + g.np.array([0, 0.95, 0.95])).all()
+        assert g.np.isclose(sliced.bounds[1], mesh.bounds[1]).all()
+        assert len(sliced.faces) == 11
+
+        # Try with more complicated mesh and make sure we get correct projections and some faces
+        origins = [bunny.bounds.mean(axis=0), bunny.bounds.mean(axis=0)+0.01*g.trimesh.unitize([1, 1, 2])]
+        normals = [g.trimesh.unitize([1, 1, 2]), -g.trimesh.unitize([1, 1, 2])]
+
+        sliced = bunny.slice_planes(plane_origins=origins,
+                                    plane_normals=normals)
+        assert len(sliced.faces) > 0
+
+        for o,n in zip(origins,normals):
+            # check the projections manually
+            dot = g.np.dot(n, (sliced.vertices - o).T)
+            # should be lots of stuff at the plane and nothing behind
+            assert g.np.isclose(dot.min(), 0.0)
 
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()