Repository: mikedh/trimesh
Pull request #801 (2020-04-18 22:02:39)

Grid Creation And GLTF Option
- Fixes bug in `trimesh.path.creation.grid`
- Breaks out tests for `trimesh.path.creation` into separate file
- Adds `merge_primitives` option to GLTF loader, which if specified combines all primitives of a GLTF `mesh` into a single `Trimesh` object (vs previously where it was one `Trimesh` object per  `primitive`). Default behavior does not change.
- reverts `Scene.add_geometry` behavior to address #802 
- releases capping of plane-sliced meshes #800 
- releases getting voxel bounds from sparse indices #795 


Patch:
diff --git a/trimesh/base.py b/trimesh/base.py
index 7a6aaae8..f1c54221 100644
--- a/trimesh/base.py
+++ b/trimesh/base.py
@@ -1982,7 +1982,8 @@ class Trimesh(Geometry):
 
     def section(self,
                 plane_normal,
-                plane_origin):
+                plane_origin,
+                **kwargs):
         """
         Returns a 3D cross section of the current mesh and a plane
         defined by origin and normal.
@@ -2007,7 +2008,8 @@ class Trimesh(Geometry):
             mesh=self,
             plane_normal=plane_normal,
             plane_origin=plane_origin,
-            return_faces=True)
+            return_faces=True,
+            **kwargs)
 
         # if the section didn't hit the mesh return None
         if len(lines) == 0:
diff --git a/trimesh/exchange/gltf.py b/trimesh/exchange/gltf.py
index 4660b904..bed54aa3 100644
--- a/trimesh/exchange/gltf.py
+++ b/trimesh/exchange/gltf.py
@@ -206,6 +206,7 @@ def export_glb(scene, extras=None, include_normals=None):
 
 def load_gltf(file_obj=None,
               resolver=None,
+              merge_primitives=False,
               **mesh_kwargs):
     """
     Load a GLTF file, which consists of a directory structure
@@ -217,6 +218,8 @@ def load_gltf(file_obj=None,
       Object containing header JSON, or None
     resolver : trimesh.visual.Resolver
       Object which can be used to load other files by name
+    merge_primitives : bool
+      If True, each GLTF 'mesh' will correspond to a single Trimesh object
     **mesh_kwargs : dict
       Passed to mesh constructor
 
@@ -245,12 +248,16 @@ def load_gltf(file_obj=None,
     # that can be used to instantiate a trimesh.Scene object
     kwargs = _read_buffers(header=tree,
                            buffers=buffers,
+                           merge_primitives=merge_primitives,
                            mesh_kwargs=mesh_kwargs,
                            resolver=resolver)
     return kwargs
 
 
-def load_glb(file_obj, resolver=None, **mesh_kwargs):
+def load_glb(file_obj,
+             resolver=None,
+             merge_primitives=False,
+             **mesh_kwargs):
     """
     Load a GLTF file in the binary GLB format into a trimesh.Scene.
 
@@ -261,6 +268,10 @@ def load_glb(file_obj, resolver=None, **mesh_kwargs):
     ------------
     file_obj : file- like object
       Containing GLB data
+    resolver : trimesh.visual.Resolver
+      Object which can be used to load other files by name
+    merge_primitives : bool
+      If True, each GLTF 'mesh' will correspond to a single Trimesh object
 
     Returns
     ------------
@@ -324,6 +335,7 @@ def load_glb(file_obj, resolver=None, **mesh_kwargs):
     # that can be used to instantiate a trimesh.Scene object
     kwargs = _read_buffers(header=header,
                            buffers=buffers,
+                           merge_primitives=merge_primitives,
                            mesh_kwargs=mesh_kwargs)
     return kwargs
 
@@ -864,7 +876,7 @@ def _parse_materials(header, views, resolver=None):
     return materials
 
 
-def _read_buffers(header, buffers, mesh_kwargs, resolver=None):
+def _read_buffers(header, buffers, mesh_kwargs, merge_primitives=False, resolver=None):
     """
     Given a list of binary data and a layout, return the
     kwargs to create a scene object.
@@ -992,9 +1004,9 @@ def _read_buffers(header, buffers, mesh_kwargs, resolver=None):
             else:
                 name = "GLTF_geometry"
 
-                # make name unique across multiple meshes
-                if len(header["meshes"]) > 1:
-                    name += "_{}".format(index)
+            # make name unique across multiple meshes
+            if name in meshes:
+                name += "_{}".format(util.unique_id())
 
             # each primitive gets it's own Trimesh object
             if len(m["primitives"]) > 1:
@@ -1006,12 +1018,44 @@ def _read_buffers(header, buffers, mesh_kwargs, resolver=None):
                 for attr in custom_attrs:
                     vertex_attributes[attr] = access[p["attributes"][attr]]
                 kwargs["vertex_attributes"] = vertex_attributes
-
             kwargs["process"] = False
-
             meshes[name] = kwargs
             mesh_prim[index].append(name)
 
+    # sometimes GLTF "meshes" come with multiple "primitives"
+    # by default we return one Trimesh object per "primitive"
+    # but if merge_primitives is True we combine the primitives
+    # for the "mesh" into a single Trimesh object
+    if merge_primitives:
+        # if we are only returning one Trimesh object
+        # replace `mesh_prim` with updated values
+        mesh_prim_replace = dict()
+        mesh_pop = []
+        for mesh_index, names in mesh_prim.items():
+            if len(names) <= 1:
+                mesh_prim_replace[mesh_index] = names
+                continue
+            name = '_'.join(names[0].split('_')[:-1])
+            mesh_pop.extend(set(names).difference({name}))
+            current = [meshes[n] for n in names]
+            v_seq = [p['vertices'] for p in current]
+            f_seq = [p['faces'] for p in current]
+            v, f = util.append_faces(v_seq, f_seq)
+            if 'metadata' in meshes[names[0]]:
+                metadata = meshes[names[0]]['metadata']
+            else:
+                metadata = {}
+            meshes[name] = {
+                'vertices': v,
+                'faces': f,
+                'metadata': metadata,
+                'process': False}
+            mesh_prim_replace[mesh_index] = [name]
+        # avoid altering inside loop
+        mesh_prim = mesh_prim_replace
+        # remove outdated meshes
+        [meshes.pop(p, None) for p in mesh_pop]
+
     # make it easier to reference nodes
     nodes = header["nodes"]
     # nodes are referenced by index
diff --git a/trimesh/intersections.py b/trimesh/intersections.py
index 04bbd237..6d197368 100644
--- a/trimesh/intersections.py
+++ b/trimesh/intersections.py
@@ -157,8 +157,8 @@ def mesh_plane(mesh,
         # 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]
+        dots = np.einsum('i,ij->j', plane_normal,
+                         (mesh.vertices - plane_origin).T)[mesh.faces]
 
     # sign of the dot product is -1, 0, or 1
     # shape is the same as mesh.faces (n,3)
@@ -421,11 +421,12 @@ def slice_faces_plane(vertices,
         Faces of source mesh to slice
     plane_normal : (3,) float
         Normal vector of plane to intersect with mesh
-    plane_origin:  (3,) float
+    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_vertices : (n, 3) float
@@ -492,22 +493,17 @@ def slice_faces_plane(vertices,
                      np.zeros((0, 3), dtype=np.int64))
             return empty
 
-        try:
-            # count the number of occurrences of each value
-            counts = np.bincount(new_faces.flatten(), minlength=len(vertices))
-            unique_verts = counts > 0
-            unique_index = np.where(unique_verts)[0]
-        except TypeError:
-            # casting failed on 32 bit windows
-            log.error('casting failed!', exc_info=True)
-            # fall back to numpy unique
-            unique_index = np.unique(new_faces.flatten())
-            # generate a mask for cumsum
-            unique_verts = np.zeros(len(vertices), dtype=np.bool)
-            unique_verts[unique_index] = True
-
-        unique_faces = (np.cumsum(unique_verts) - 1)[new_faces]
-        return vertices[unique_index], unique_faces
+        # find the unique indices in the new faces
+        # using an integer-only unique function
+        unique, inverse = grouping.unique_bincount(new_faces.reshape(-1),
+                                                   minlength=len(vertices),
+                                                   return_inverse=True)
+
+        # use the unique indices for our final vertices and faces
+        final_vert = vertices[unique]
+        final_face = inverse.reshape((-1, 3))
+
+        return final_vert, final_face
 
     # Extract the intersections of each triangle's edges with the plane
     o = cut_triangles                               # origins
@@ -586,7 +582,7 @@ def slice_faces_plane(vertices,
         new_faces = np.append(new_faces, new_tri_faces, axis=0)
 
     # find the unique indices in the new faces
-    # using an integer- only unique function
+    # using an integer-only unique function
     unique, inverse = grouping.unique_bincount(new_faces.reshape(-1),
                                                minlength=len(new_vertices),
                                                return_inverse=True)
@@ -601,6 +597,8 @@ def slice_faces_plane(vertices,
 def slice_mesh_plane(mesh,
                      plane_normal,
                      plane_origin,
+                     cap=False,
+                     cached_dots=None,
                      **kwargs):
     """
     Slice a mesh with a plane, returning a new mesh that is the
@@ -612,11 +610,10 @@ def slice_mesh_plane(mesh,
       Source mesh to slice
     plane_normal : (3,) float
       Normal vector of plane to intersect with mesh
-    plane_origin:  (3,) float
+    plane_origin :  (3,) float
       Point on plane to intersect with mesh
-    cap: bool
+    cap : bool
       If True, cap the result with a triangulated polygon
-
     cached_dots : (n, 3) float
         If an external function has stored dot
         products pass them here to avoid recomputing
@@ -632,6 +629,7 @@ def slice_mesh_plane(mesh,
 
     # avoid circular import
     from .base import Trimesh
+    from .creation import triangulate_polygon
 
     # check input plane
     plane_normal = np.asanyarray(plane_normal,
@@ -648,22 +646,67 @@ def slice_mesh_plane(mesh,
     if not shape_ok:
         raise ValueError('plane origins and normals must be (n, 3)!')
 
-    # start with original vertices and faces
+    # start with copy of original mesh, faces, and vertices
+    sliced_mesh = mesh.copy()
     vertices = mesh.vertices.copy()
     faces = mesh.faces.copy()
 
     # slice away specified planes
     for origin, normal in zip(plane_origin.reshape((-1, 3)),
                               plane_normal.reshape((-1, 3))):
+
+        # calculate dots here if not passed in to save time
+        # in case of cap
+        if cached_dots is None:
+            # 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 faces (n,3)
+            dots = np.einsum('i,ij->j', normal,
+                             (vertices - origin).T)[faces]
+        else:
+            dots = cached_dots
         # save the new vertices and faces
         vertices, faces = slice_faces_plane(vertices=vertices,
                                             faces=faces,
                                             plane_normal=normal,
                                             plane_origin=origin,
-                                            **kwargs)
-
-    # create a mesh from the sliced result
-    new_mesh = Trimesh(vertices=vertices,
-                       faces=faces,
-                       process=False)
-    return new_mesh
+                                            cached_dots=dots)
+
+        # check if cap arg specified
+        if cap:
+            # check if mesh is watertight (can't cap if not)
+            if not sliced_mesh.is_watertight:
+                raise ValueError('Input mesh must be watertight to cap slice')
+
+            path = sliced_mesh.section(plane_normal=normal,
+                                       plane_origin=origin,
+                                       cached_dots=dots)
+
+            # transform Path3D onto XY plane for triangulation
+            on_plane, to_3D = path.to_planar()
+
+            # triangulate each closed region of 2D cap
+            # without adding any new vertices
+            v, f = [], []
+            for polygon in on_plane.polygons_full:
+                t = triangulate_polygon(
+                    polygon, triangle_args='p', allow_boundary_steiner=False)
+                v.append(t[0])
+                f.append(t[1])
+
+            # append regions and reindex
+            vf, ff = util.append_faces(v, f)
+
+            # make vertices 3D and transform back to mesh frame
+            vf = np.column_stack((vf, np.zeros(len(vf))))
+            vf = transformations.transform_points(vf, to_3D)
+
+            # add cap vertices and faces and reindex
+            vertices, faces = util.append_faces([vertices, vf], [faces, ff])
+
+            # Update mesh with cap (processing needed to merge vertices)
+            sliced_mesh = Trimesh(vertices=vertices, faces=faces)
+            vertices, faces = sliced_mesh.vertices.copy(), sliced_mesh.faces.copy()
+
+    # return the sliced mesh
+    return Trimesh(vertices=vertices, faces=faces, process=False)
diff --git a/trimesh/path/creation.py b/trimesh/path/creation.py
index 6bcd8a66..d24bd0b2 100644
--- a/trimesh/path/creation.py
+++ b/trimesh/path/creation.py
@@ -2,6 +2,7 @@ import numpy as np
 
 from . import arc
 from .entities import Line, Arc
+from ..geometry import plane_transform
 
 from .. import util
 from .. import transformations
@@ -300,7 +301,7 @@ def grid(side,
 
     # if plane was passed instead of transform create the matrix here
     if (transform is None and plane_origin is not None and plane_normal is not None):
-        transform = np.linalg.inv(transformations.plane_transform(
+        transform = np.linalg.inv(plane_transform(
             origin=plane_origin, normal=plane_normal))
 
     # stack vertices to 3D
diff --git a/trimesh/primitives.py b/trimesh/primitives.py
index c3adeca7..799c53af 100644
--- a/trimesh/primitives.py
+++ b/trimesh/primitives.py
@@ -612,7 +612,7 @@ class Box(_Primitive):
         outline : trimesh.path.Path3D
           Outline of box primitive
         """
-        # do the import in function to keep soft dependancy
+        # do the import in function to keep soft dependency
         from .path.creation import box_outline
         # return outline with same size as primitive
         return box_outline(
diff --git a/trimesh/scene/scene.py b/trimesh/scene/scene.py
index 5c88b56c..a19c0181 100644
--- a/trimesh/scene/scene.py
+++ b/trimesh/scene/scene.py
@@ -127,21 +127,12 @@ class Scene(Geometry):
         # PointCloud objects will look like a sequence
         elif util.is_sequence(geometry):
             # if passed a sequence add all elements
-            for i, value in enumerate(geometry):
-                if i == 0:
-                    node_name = self.add_geometry(
-                        geometry=value,
-                        node_name=node_name,
-                        geom_name=geom_name,
-                        parent_node_name=parent_node_name,
-                        transform=transform)
-                else:
-                    self.add_geometry(
-                        geometry=value,
-                        geom_name=geom_name,
-                        parent_node_name=node_name)
-            return
-
+            return [self.add_geometry(
+                geometry=value,
+                node_name=node_name,
+                geom_name=geom_name,
+                parent_node_name=parent_node_name,
+                transform=transform) for value in geometry]
         elif isinstance(geometry, dict):
             # if someone passed us a dict of geometry
             for key, value in geometry.items():
diff --git a/trimesh/version.py b/trimesh/version.py
index 8b74917b..d092c3d2 100644
--- a/trimesh/version.py
+++ b/trimesh/version.py
@@ -1,1 +1,1 @@
-__version__ = '3.6.30'
+__version__ = '3.6.32'
diff --git a/trimesh/voxel/base.py b/trimesh/voxel/base.py
index a517979a..f137c33f 100644
--- a/trimesh/voxel/base.py
+++ b/trimesh/voxel/base.py
@@ -11,6 +11,7 @@ from . import transforms
 from . import morphology
 from . import encoding as enc
 
+from .. import bounds as bounds_module
 from .. import caching
 from .. import transformations as tr
 
@@ -141,8 +142,14 @@ class VoxelGrid(Geometry):
 
     @caching.cache_decorator
     def bounds(self):
-        points = self.points
-        bounds = np.array([points.min(axis=0), points.max(axis=0)])
+        indices = self.sparse_indices
+        # get all 8 corners of the AABB
+        corners = bounds_module.corners([indices.min(axis=0) - 0.5,
+                                         indices.max(axis=0) + 0.5])
+        # transform these corners to a new frame
+        corners = self._transform.transform_points(corners)
+        # get the AABB of corners in-frame
+        bounds = np.array([corners.min(axis=0), corners.max(axis=0)])
         bounds.flags.writeable = False
         return bounds
 


Test patch:
diff --git a/tests/test_binvox.py b/tests/test_binvox.py
index cd8c0799..9c1077d1 100644
--- a/tests/test_binvox.py
+++ b/tests/test_binvox.py
@@ -29,7 +29,9 @@ class BinvoxTest(g.unittest.TestCase):
             [0, 0, 0, 1]
         ]))
         dense = dense.transpose((0, 2, 1))
-        np.testing.assert_allclose(base.bounds, [translate, translate + scale])
+        bound_min = translate - 0.5 * s
+        bound_max = translate + scale + 0.5 * s
+        np.testing.assert_allclose(base.bounds, [bound_min, bound_max])
         np.testing.assert_equal(base.encoding.dense, dense)
 
         if binvox.binvox_encoder is None:
diff --git a/tests/test_gltf.py b/tests/test_gltf.py
index e39f25d7..00bfda7a 100644
--- a/tests/test_gltf.py
+++ b/tests/test_gltf.py
@@ -103,7 +103,7 @@ class GLTFTest(g.unittest.TestCase):
         m = g.get_mesh('testplate.glb')
         assert m.units == 'meters'
 
-    def test_gltf(self):
+    def test_basic(self):
         # split a multibody mesh into a scene
         scene = g.trimesh.scene.split_scene(
             g.get_mesh('cycloidal.ply'))
@@ -137,7 +137,7 @@ class GLTFTest(g.unittest.TestCase):
             # will assert round trip is roughly equal
             g.scene_equal(rd, scene)
 
-    def test_gltf_merge_buffers(self):
+    def test_merge_buffers(self):
         # split a multibody mesh into a scene
         scene = g.trimesh.scene.split_scene(
             g.get_mesh('cycloidal.ply'))
@@ -157,7 +157,16 @@ class GLTFTest(g.unittest.TestCase):
         # check to make sure the geometry keys are the same
         assert set(reloaded.geometry.keys()) == set(scene.geometry.keys())
 
-    def test_gltf_optional_camera(self):
+    def test_merge_primitives(self):
+        # test to see if the `merge_primitives` logic is working
+        a = g.get_mesh('CesiumMilkTruck.glb')
+        assert len(a.geometry) == 4
+
+        # should combine the multiple primitives into a single mesh
+        b = g.get_mesh('CesiumMilkTruck.glb', merge_primitives=True)
+        assert len(b.geometry) == 2
+
+    def test_optional_camera(self):
         gltf_cameras_key = 'cameras'
 
         # if there's no camera in the scene, then it shouldn't be added to the gltf
diff --git a/tests/test_path_creation.py b/tests/test_path_creation.py
new file mode 100644
index 00000000..3389f22b
--- /dev/null
+++ b/tests/test_path_creation.py
@@ -0,0 +1,75 @@
+try:
+    from . import generic as g
+except BaseException:
+    import generic as g
+
+
+class CreationTests(g.unittest.TestCase):
+
+    def test_circle_pattern(self):
+        from trimesh.path import creation
+        pattern = creation.circle_pattern(pattern_radius=1.0,
+                                          circle_radius=0.1,
+                                          count=4)
+        assert len(pattern.entities) == 4
+        assert len(pattern.polygons_closed) == 4
+        assert len(pattern.polygons_full) == 4
+
+        # should be a valid Path2D
+        g.check_path2D(pattern)
+
+    def test_circle(self):
+        from trimesh.path import creation
+        circle = creation.circle(radius=1.0, center=(1.0, 1.0))
+
+        # it's a discrete circle
+        assert g.np.isclose(circle.area, g.np.pi, rtol=0.01)
+        # should be centered at 0
+        assert g.np.allclose(
+            circle.polygons_full[0].centroid, [
+                1.0, 1.0], atol=1e-3)
+        assert len(circle.entities) == 1
+        assert len(circle.polygons_closed) == 1
+        assert len(circle.polygons_full) == 1
+
+        # should be a valid Path2D
+        g.check_path2D(circle)
+
+    def test_rect(self):
+        from trimesh.path import creation
+
+        # create a single rectangle
+        pattern = creation.rectangle([[0, 0], [2, 3]])
+        assert len(pattern.entities) == 1
+        assert len(pattern.polygons_closed) == 1
+        assert len(pattern.polygons_full) == 1
+        assert g.np.isclose(pattern.area, 6.0)
+        # should be a valid Path2D
+        g.check_path2D(pattern)
+
+        # make 10 untouching rectangles
+        pattern = creation.rectangle(
+            g.np.arange(40).reshape((-1, 2, 2)))
+        assert len(pattern.entities) == 10
+        assert len(pattern.polygons_closed) == 10
+        assert len(pattern.polygons_full) == 10
+        # should be a valid Path2D
+        g.check_path2D(pattern)
+
+    def test_grid(self):
+        grid = g.trimesh.path.creation.grid(side=5.0)
+        print(grid.extents)
+        assert g.np.allclose(grid.extents, [10, 10, 0])
+
+        # check grid along a plane
+        grid = g.trimesh.path.creation.grid(side=10.0,
+                                            plane_origin=[5.0, 0, 0],
+                                            plane_normal=[1, 0, 0])
+        # make sure plane is applied correctly
+        assert g.np.allclose(grid.extents, [0, 20, 20])
+        assert g.np.allclose(grid.bounds, [[5, -10, -10], [5, 10, 10]])
+
+
+if __name__ == '__main__':
+    g.trimesh.util.attach_to_log()
+    g.unittest.main()
diff --git a/tests/test_paths.py b/tests/test_paths.py
index b20d65cd..27a39252 100644
--- a/tests/test_paths.py
+++ b/tests/test_paths.py
@@ -338,59 +338,6 @@ class SectionTest(g.unittest.TestCase):
             g.check_path2D(planar)
 
 
-class CreationTests(g.unittest.TestCase):
-
-    def test_circle_pattern(self):
-        from trimesh.path import creation
-        pattern = creation.circle_pattern(pattern_radius=1.0,
-                                          circle_radius=0.1,
-                                          count=4)
-        assert len(pattern.entities) == 4
-        assert len(pattern.polygons_closed) == 4
-        assert len(pattern.polygons_full) == 4
-
-        # should be a valid Path2D
-        g.check_path2D(pattern)
-
-    def test_circle(self):
-        from trimesh.path import creation
-        circle = creation.circle(radius=1.0, center=(1.0, 1.0))
-
-        # it's a discrete circle
-        assert g.np.isclose(circle.area, g.np.pi, rtol=0.01)
-        # should be centered at 0
-        assert g.np.allclose(
-            circle.polygons_full[0].centroid, [
-                1.0, 1.0], atol=1e-3)
-        assert len(circle.entities) == 1
-        assert len(circle.polygons_closed) == 1
-        assert len(circle.polygons_full) == 1
-
-        # should be a valid Path2D
-        g.check_path2D(circle)
-
-    def test_rect(self):
-        from trimesh.path import creation
-
-        # create a single rectangle
-        pattern = creation.rectangle([[0, 0], [2, 3]])
-        assert len(pattern.entities) == 1
-        assert len(pattern.polygons_closed) == 1
-        assert len(pattern.polygons_full) == 1
-        assert g.np.isclose(pattern.area, 6.0)
-        # should be a valid Path2D
-        g.check_path2D(pattern)
-
-        # make 10 untouching rectangles
-        pattern = creation.rectangle(
-            g.np.arange(40).reshape((-1, 2, 2)))
-        assert len(pattern.entities) == 10
-        assert len(pattern.polygons_closed) == 10
-        assert len(pattern.polygons_full) == 10
-        # should be a valid Path2D
-        g.check_path2D(pattern)
-
-
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()
     g.unittest.main()
diff --git a/tests/test_scene.py b/tests/test_scene.py
index 0cb72103..ea484260 100644
--- a/tests/test_scene.py
+++ b/tests/test_scene.py
@@ -170,7 +170,7 @@ class SceneTests(g.unittest.TestCase):
         #     scene.add_geometry(geometry[1], parent_node_name='voxels')
         scene = g.trimesh.Scene()
         geometry = [g.trimesh.creation.box(), g.trimesh.creation.box()]
-        scene.add_geometry(geometry, node_name='voxels')
+        scene.add_geometry(geometry)
         assert len(scene.graph.nodes_geometry) == 2
 
     def test_add_concat(self):
diff --git a/tests/test_section.py b/tests/test_section.py
index ab72fd8f..646016da 100644
--- a/tests/test_section.py
+++ b/tests/test_section.py
@@ -140,7 +140,7 @@ class SliceTest(g.unittest.TestCase):
         mesh = g.trimesh.creation.box()
 
         # Cut corner off of box and make sure the bounds and number of faces is correct
-        # (Tests new triangles, but not new quads or triangles contained entirely)
+        # Tests new triangles, but not new quads or triangles contained entirely
         plane_origin = mesh.bounds[1] - 0.05
         plane_normal = mesh.bounds[1]
 
@@ -208,6 +208,94 @@ class SliceTest(g.unittest.TestCase):
             # should be lots of stuff at the plane and nothing behind
             assert g.np.isclose(dot.min(), 0.0)
 
+        # Test cap on more complicated watertight mesh to make sure the
+        # resulting mesh is still watertight and slice is correct
+        featuretype = g.get_mesh('featuretype.STL')
+
+        origins = [featuretype.center_mass, featuretype.center_mass
+                   + 0.01 * g.trimesh.unitize([1, 0, 2])]
+        normals = [g.trimesh.unitize([1, 1, 1]), g.trimesh.unitize([1, 2, 3])]
+
+    def test_cap(self):
+
+        try:
+            from triangle import triangulate  # NOQA
+        except BaseException as E:
+            if g.all_dep:
+                raise E
+            else:
+                return
+
+        mesh = g.trimesh.creation.box()
+
+        # Cut corner off of box and make sure the bounds and number of faces is correct
+        # Tests new triangles, but not new quads or triangles contained entirely
+        plane_origin = mesh.bounds[1] - 0.05
+        plane_normal = mesh.bounds[1]
+
+        # Same test with capping (should only add three more triangles)
+        sliced_capped = mesh.slice_plane(plane_origin=plane_origin,
+                                         plane_normal=plane_normal,
+                                         cap=True)
+
+        assert len(sliced_capped.faces) == 8
+        assert g.np.isclose(sliced_capped.bounds[0], mesh.bounds[1] - 0.15).all()
+        assert g.np.isclose(sliced_capped.bounds[1], mesh.bounds[1]).all()
+        assert sliced_capped.is_watertight
+
+        # Cut top off of box and make sure bounds and number of faces is correct
+        # Tests new quads and entirely contained triangles
+        plane_origin = mesh.bounds[1] - 0.05
+        plane_normal = g.np.array([0, 0, 1])
+
+        # Same test with capping (should only add six triangles)
+        sliced_capped = mesh.slice_plane(plane_origin=plane_origin,
+                                         plane_normal=plane_normal,
+                                         cap=True)
+
+        assert len(sliced_capped.faces) == 20
+        assert g.np.isclose(
+            sliced_capped.bounds[0], mesh.bounds[0] + g.np.array([0, 0, 0.95])).all()
+        assert g.np.isclose(sliced_capped.bounds[1], mesh.bounds[1]).all()
+        assert sliced_capped.is_watertight
+
+        # 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])]
+
+        # Test cap for multiple slices to check watertightness
+        # (should add nine triangles)
+        sliced_capped = mesh.slice_plane(plane_origin=plane_origins,
+                                         plane_normal=plane_normals,
+                                         cap=True)
+
+        assert len(sliced_capped.faces) == 20
+        assert g.np.isclose(
+            sliced_capped.bounds[0], mesh.bounds[0] + g.np.array([0, 0.95, 0.95])).all()
+        assert g.np.isclose(sliced_capped.bounds[1], mesh.bounds[1]).all()
+        assert sliced_capped.is_watertight
+
+        # Test cap on more complicated watertight mesh to make sure the
+        # resulting mesh is still watertight and slice is correct
+        featuretype = g.get_mesh('featuretype.STL')
+
+        origins = [featuretype.center_mass, featuretype.center_mass
+                   + 0.01 * g.trimesh.unitize([1, 0, 2])]
+        normals = [g.trimesh.unitize([1, 1, 1]), g.trimesh.unitize([1, 2, 3])]
+
+        sliced_capped = featuretype.slice_plane(plane_origin=origins,
+                                                plane_normal=normals,
+                                                cap=True)
+
+        assert len(sliced_capped.faces) > 0
+        assert sliced_capped.is_watertight
+        for o, n in zip(origins, normals):
+            # check the projections manually
+            dot = g.np.dot(n, (sliced_capped.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()