Repository: mikedh/trimesh
Pull request #2319 (2024-11-05 20:10:50)

Release: Iteration Tests 
- make `trimesh.iteration` standard library only since it was only using `numpy.log2`
- separate `itertools`-like methods tests into `test_iteration.py`
- release #2320 
- release #2322 
- release #2313 


Patch:
diff --git a/pyproject.toml b/pyproject.toml
index b385a492..3609dcd9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,7 +5,7 @@ requires = ["setuptools >= 61.0", "wheel"]
 [project]
 name = "trimesh"
 requires-python = ">=3.8"
-version = "4.5.2"
+version = "4.5.3"
 authors = [{name = "Michael Dawson-Haggerty", email = "mikedh@kerfed.com"}]
 license = {file = "LICENSE.md"}
 description = "Import, export, process, analyze and view triangular meshes."
diff --git a/trimesh/__init__.py b/trimesh/__init__.py
index 0742d68f..40d3b85c 100644
--- a/trimesh/__init__.py
+++ b/trimesh/__init__.py
@@ -82,9 +82,8 @@ except BaseException as E:
 
 
 __all__ = [
-    "PointCloud",
     "Geometry",
-    "Trimesh",
+    "PointCloud",
     "Scene",
     "Trimesh",
     "__version__",
@@ -103,8 +102,8 @@ __all__ = [
     "graph",
     "grouping",
     "inertia",
-    "iteration",
     "intersections",
+    "iteration",
     "load",
     "load_mesh",
     "load_path",
diff --git a/trimesh/exchange/binvox.py b/trimesh/exchange/binvox.py
index bc6e10aa..67eadd5a 100644
--- a/trimesh/exchange/binvox.py
+++ b/trimesh/exchange/binvox.py
@@ -201,7 +201,9 @@ def voxel_from_binvox(rle_data, shape, translate=None, scale=1.0, axis_order="xz
     elif axis_order is None or axis_order == "xyz":
         encoding = encoding.reshape(shape)
     else:
-        raise ValueError("Invalid axis_order '%s': must be None, 'xyz' or 'xzy'")
+        raise ValueError(
+            "Invalid axis_order '%s': must be None, 'xyz' or 'xzy'", axis_order
+        )
 
     assert encoding.shape == shape
 
@@ -423,7 +425,7 @@ class Binvoxer:
             )
 
         if dimension > 1024 and not exact:
-            raise ValueError("Maximum dimension using exact is 1024, got %d" % dimension)
+            raise ValueError("Maximum dimension using exact is 1024, got %d", dimension)
         if file_type not in Binvoxer.SUPPORTED_OUTPUT_TYPES:
             raise ValueError(
                 f"file_type {file_type} not in set of supported output types {Binvoxer.SUPPORTED_OUTPUT_TYPES!s}"
@@ -471,7 +473,7 @@ class Binvoxer:
             times = np.log2(downsample_factor)
             if int(times) != times:
                 raise ValueError(
-                    "downsample_factor must be a power of 2, got %d" % downsample_factor
+                    "downsample_factor must be a power of 2, got %d", downsample_factor
                 )
             args.extend(("-down",) * int(times))
         if downsample_threshold is not None:
diff --git a/trimesh/exchange/gltf.py b/trimesh/exchange/gltf.py
index 89ab3b0c..63217877 100644
--- a/trimesh/exchange/gltf.py
+++ b/trimesh/exchange/gltf.py
@@ -1029,14 +1029,14 @@ def _build_accessor(array):
         if vec_length > 4:
             raise ValueError("The GLTF spec does not support vectors larger than 4")
         if vec_length > 1:
-            data_type = "VEC%d" % vec_length
+            data_type = f"VEC{int(vec_length)}"
         else:
             data_type = "SCALAR"
 
     if len(shape) == 3:
         if shape[2] not in [2, 3, 4]:
             raise ValueError("Matrix types must have 4, 9 or 16 components")
-        data_type = "MAT%d" % shape[2]
+        data_type = f"MAT{int(shape[2])}"
 
     # get the array data type as a str stripping off endian
     lookup = array.dtype.str.lstrip("<>")
diff --git a/trimesh/interfaces/__init__.py b/trimesh/interfaces/__init__.py
index 45416ec8..4f58ef4a 100644
--- a/trimesh/interfaces/__init__.py
+++ b/trimesh/interfaces/__init__.py
@@ -1,4 +1,4 @@
 from . import blender, gmsh
 
 # add to __all__ as per pep8
-__all__ = ["gmsh", "blender"]
+__all__ = ["blender", "gmsh"]
diff --git a/trimesh/iteration.py b/trimesh/iteration.py
index 5994341f..ca1e7080 100644
--- a/trimesh/iteration.py
+++ b/trimesh/iteration.py
@@ -1,4 +1,4 @@
-import numpy as np
+from math import log2
 
 from .typed import Any, Callable, Iterable, List, NDArray, Sequence, Union
 
@@ -54,12 +54,17 @@ def reduce_cascade(operation: Callable, items: Union[Sequence, NDArray]):
         # skip the loop overhead for a single pair
         return operation(items[0], items[1])
 
-    for _ in range(int(1 + np.log2(len(items)))):
+    for _ in range(int(1 + log2(len(items)))):
         results = []
-        for i in np.arange(len(items) // 2) * 2:
+
+        # loop over pairs of items.
+        items_mod = len(items) % 2
+        for i in range(0, len(items) - items_mod, 2):
             results.append(operation(items[i], items[i + 1]))
 
-        if len(items) % 2:
+        # if we had a non-even number of items it will have been
+        # skipped by the loop so append it to our list
+        if items_mod != 0:
             results.append(items[-1])
 
         items = results
@@ -117,7 +122,7 @@ def chain(*args: Union[Iterable[Any], Any, None]) -> List[Any]:
     # extend if it's a sequence, otherwise append
     [
         chained.extend(a)
-        if (hasattr(a, "__iter__") and not isinstance(a, str))
+        if (hasattr(a, "__iter__") and not isinstance(a, (str, bytes)))
         else chained.append(a)
         for a in args
         if a is not None
diff --git a/trimesh/path/exchange/load.py b/trimesh/path/exchange/load.py
index b6b28fc3..16dea51e 100644
--- a/trimesh/path/exchange/load.py
+++ b/trimesh/path/exchange/load.py
@@ -1,50 +1,55 @@
-import os
-
 from ... import util
 from ...exchange.ply import load_ply
+from ...typed import Optional
 from ..path import Path
 from . import misc
 from .dxf import _dxf_loaders
 from .svg_io import svg_to_path
 
 
-def load_path(file_obj, file_type=None, **kwargs):
+def load_path(file_obj, file_type: Optional[str] = None, **kwargs):
     """
     Load a file to a Path file_object.
 
     Parameters
     -----------
-    file_obj : One of the following:
+    file_obj
+      Accepts many types:
          - Path, Path2D, or Path3D file_objects
          - open file file_object (dxf or svg)
          - file name (dxf or svg)
          - shapely.geometry.Polygon
          - shapely.geometry.MultiLineString
          - dict with kwargs for Path constructor
-         - (n,2,(2|3)) float, line segments
-    file_type : str
+         - `(n, 2, (2|3)) float` line segments
+    file_type
         Type of file is required if file
-        file_object passed.
+        object is passed.
 
     Returns
     ---------
     path : Path, Path2D, Path3D file_object
-        Data as a native trimesh Path file_object
+      Data as a native trimesh Path file_object
     """
     # avoid a circular import
     from ...exchange.load import load_kwargs
 
+    if isinstance(file_type, str):
+        # we accept full file names here so make sure we
+        file_type = util.split_extension(file_type).lower()
+
     # record how long we took
     tic = util.now()
 
     if isinstance(file_obj, Path):
-        # we have been passed a Path file_object so
-        # do nothing and return the passed file_object
+        # we have been passed a file object that is already a loaded
+        # trimesh.path.Path object so do nothing and return
         return file_obj
     elif util.is_file(file_obj):
         # for open file file_objects use loaders
         if file_type == "ply":
-            # we cannot register this exporter to path_loaders since this is already reserved by TriMesh in ply format in trimesh.load()
+            # we cannot register this exporter to path_loaders
+            # since this is already reserved for 3D values in `trimesh.load`
             kwargs.update(load_ply(file_obj, file_type=file_type))
         else:
             kwargs.update(path_loaders[file_type](file_obj, file_type=file_type))
@@ -52,7 +57,7 @@ def load_path(file_obj, file_type=None, **kwargs):
         # strings passed are evaluated as file file_objects
         with open(file_obj, "rb") as f:
             # get the file type from the extension
-            file_type = os.path.splitext(file_obj)[-1][1:].lower()
+            file_type = util.split_extension(file_obj).lower()
             if file_type == "ply":
                 # we cannot register this exporter to path_loaders since this is already reserved by TriMesh in ply format in trimesh.load()
                 kwargs.update(load_ply(f, file_type=file_type))
diff --git a/trimesh/path/packing.py b/trimesh/path/packing.py
index 8af41638..2e0ff619 100644
--- a/trimesh/path/packing.py
+++ b/trimesh/path/packing.py
@@ -8,7 +8,7 @@ Pack rectangular regions onto larger rectangular regions.
 import numpy as np
 
 from ..constants import log, tol
-from ..typed import Integer, Number, Optional
+from ..typed import ArrayLike, Integer, NDArray, Number, Optional, float64
 from ..util import allclose, bounds_tree
 
 # floating point zero
@@ -692,7 +692,7 @@ def visualize(extents, bounds):
     return Scene(meshes)
 
 
-def roll_transform(bounds, extents):
+def roll_transform(bounds: ArrayLike, extents: ArrayLike) -> NDArray[float64]:
     """
     Packing returns rotations with integer "roll" which
     needs to be converted into a homogeneous rotation matrix.
diff --git a/trimesh/path/util.py b/trimesh/path/util.py
index e9862ae9..196f81c0 100644
--- a/trimesh/path/util.py
+++ b/trimesh/path/util.py
@@ -3,7 +3,7 @@ import numpy as np
 from ..util import is_ccw  # NOQA
 
 
-def concatenate(paths):
+def concatenate(paths, **kwargs):
     """
     Concatenate multiple paths into a single path.
 
@@ -11,6 +11,8 @@ def concatenate(paths):
     -------------
     paths : (n,) Path
       Path objects to concatenate
+    kwargs
+      Passed through to the path constructor
 
     Returns
     -------------
@@ -52,6 +54,6 @@ def concatenate(paths):
     # generate the single new concatenated path
     # use input types so we don't have circular imports
     concat = type(path)(
-        metadata=metadata, entities=entities, vertices=np.vstack(vertices)
+        metadata=metadata, entities=entities, vertices=np.vstack(vertices), **kwargs
     )
     return concat
diff --git a/trimesh/remesh.py b/trimesh/remesh.py
index 29293ff5..cbfddbcd 100644
--- a/trimesh/remesh.py
+++ b/trimesh/remesh.py
@@ -90,11 +90,7 @@ def subdivide(
     if vertex_attributes is not None:
         new_attributes = {}
         for key, values in vertex_attributes.items():
-            attr_tris = values[faces_subset]
-            attr_mid = np.vstack(
-                [attr_tris[:, g, :].mean(axis=1) for g in [[0, 1], [1, 2], [2, 0]]]
-            )
-            attr_mid = attr_mid[unique]
+            attr_mid = values[edges[unique]].mean(axis=1)
             new_attributes[key] = np.vstack((values, attr_mid))
         return new_vertices, new_faces, new_attributes
 
diff --git a/trimesh/typed.py b/trimesh/typed.py
index b27698b5..4118dfe8 100644
--- a/trimesh/typed.py
+++ b/trimesh/typed.py
@@ -20,9 +20,20 @@ if version_info >= (3, 9):
     List = list
     Tuple = tuple
     Dict = dict
+    Set = set
     from collections.abc import Callable, Hashable, Iterable, Mapping, Sequence
 else:
-    from typing import Callable, Dict, Hashable, Iterable, List, Mapping, Sequence, Tuple
+    from typing import (
+        Callable,
+        Dict,
+        Hashable,
+        Iterable,
+        List,
+        Mapping,
+        Sequence,
+        Set,
+        Tuple,
+    )
 
 # most loader routes take `file_obj` which can either be
 # a file-like object or a file path, or sometimes a dict
@@ -51,20 +62,21 @@ __all__ = [
     "Any",
     "ArrayLike",
     "BinaryIO",
+    "Callable",
     "Dict",
+    "Hashable",
     "Integer",
     "Iterable",
     "List",
     "Loadable",
+    "Mapping",
     "NDArray",
     "Number",
     "Optional",
     "Sequence",
+    "Set",
     "Stream",
     "Tuple",
     "float64",
     "int64",
-    "Mapping",
-    "Callable",
-    "Hashable",
 ]
diff --git a/trimesh/util.py b/trimesh/util.py
index 66ac4775..38fdbd7b 100644
--- a/trimesh/util.py
+++ b/trimesh/util.py
@@ -1,5 +1,5 @@
 """
-"Grab bag" of utility functions.
+Grab bag of utility functions.
 """
 
 import abc
@@ -25,7 +25,7 @@ import numpy as np
 from .iteration import chain
 
 # use our wrapped types for wider version compatibility
-from .typed import Union
+from .typed import Dict, Iterable, Optional, Set, Union
 
 # create a default logger
 log = logging.getLogger("trimesh")
@@ -1906,7 +1906,7 @@ def compress(info, **kwargs):
     return compressed
 
 
-def split_extension(file_name, special=None):
+def split_extension(file_name, special=None) -> str:
     """
     Find the file extension of a file name, including support for
     special case multipart file extensions (like .tar.gz)
@@ -2359,7 +2359,11 @@ def is_ccw(points, return_all=False):
     return ccw, area, centroid
 
 
-def unique_name(start, contains, counts=None):
+def unique_name(
+    start: Optional[str],
+    contains: Union[Set, Mapping, Iterable],
+    counts: Optional[Dict] = None,
+):
     """
     Deterministically generate a unique name not
     contained in a dict, set or other grouping with
diff --git a/trimesh/voxel/encoding.py b/trimesh/voxel/encoding.py
index 805cf34a..de315963 100644
--- a/trimesh/voxel/encoding.py
+++ b/trimesh/voxel/encoding.py
@@ -768,17 +768,20 @@ class ShapedEncoding(LazyIndexMap):
             size = np.abs(size)
             if self._data.size % size != 0:
                 raise ValueError(
-                    "cannot reshape encoding of size %d into shape %s"
-                    % (self._data.size, str(self._shape))
+                    "cannot reshape encoding of size %d into shape %s",
+                    self._data.size,
+                    str(self._shape),
                 )
+
             rem = self._data.size // size
             self._shape = tuple(rem if s == -1 else s for s in self._shape)
         elif nn > 2:
             raise ValueError("shape cannot have more than one -1 value")
         elif np.prod(self._shape) != self._data.size:
             raise ValueError(
-                "cannot reshape encoding of size %d into shape %s"
-                % (self._data.size, str(self._shape))
+                "cannot reshape encoding of size %d into shape %s",
+                self._data.size,
+                str(self._shape),
             )
 
     def _from_base_indices(self, base_indices):
@@ -818,9 +821,11 @@ class TransposedEncoding(LazyIndexMap):
             raise ValueError(f"base_encoding must be an Encoding, got {base_encoding!s}")
         if len(base_encoding.shape) != len(perm):
             raise ValueError(
-                "base_encoding has %d ndims - cannot transpose with perm %s"
-                % (base_encoding.ndims, str(perm))
+                "base_encoding has %d ndims - cannot transpose with perm %s",
+                base_encoding.ndims,
+                str(perm),
             )
+
         super().__init__(base_encoding)
         perm = np.array(perm, dtype=np.int64)
         if not all(i in perm for i in range(base_encoding.ndims)):
@@ -898,7 +903,7 @@ class FlippedEncoding(LazyIndexMap):
         super().__init__(encoding)
         if not all(0 <= a < self._data.ndims for a in axes):
             raise ValueError(
-                "Invalid axes %s for %d-d encoding" % (str(axes), self._data.ndims)
+                "Invalid axes %s for %d-d encoding", str(axes), self._data.ndims
             )
 
     def _to_base_indices(self, indices):
diff --git a/trimesh/voxel/morphology.py b/trimesh/voxel/morphology.py
index dacdf135..38692cca 100644
--- a/trimesh/voxel/morphology.py
+++ b/trimesh/voxel/morphology.py
@@ -42,7 +42,7 @@ def _sparse_indices(encoding, rank=None):
 
 def _assert_rank(value, rank):
     if len(value.shape) != rank:
-        raise ValueError("Expected rank %d, got shape %s" % (rank, str(value.shape)))
+        raise ValueError("Expected rank %d, got shape %s", rank, str(value.shape))
 
 
 def _assert_sparse_rank(value, rank=None):
@@ -51,7 +51,7 @@ def _assert_sparse_rank(value, rank=None):
     if rank is not None:
         if value.shape[-1] != rank:
             raise ValueError(
-                "sparse_indices.shape[1] must be %d, got %d" % (rank, value.shape[-1])
+                "sparse_indices.shape[1] must be %d, got %d", rank, value.shape[-1]
             )
 
 
diff --git a/trimesh/voxel/ops.py b/trimesh/voxel/ops.py
index 8e1efed8..eb910b7d 100644
--- a/trimesh/voxel/ops.py
+++ b/trimesh/voxel/ops.py
@@ -2,6 +2,7 @@ import numpy as np
 
 from .. import util
 from ..constants import log
+from ..typed import ArrayLike, Number, Optional, Union
 
 
 def fill_orthographic(dense):
@@ -95,7 +96,11 @@ def fill_base(sparse_indices):
 fill_voxelization = fill_base
 
 
-def matrix_to_marching_cubes(matrix, pitch=1.0):
+def matrix_to_marching_cubes(
+    matrix: ArrayLike,
+    pitch: Union[Number, ArrayLike] = 1.0,
+    threshold: Optional[Number] = None,
+):
     """
     Convert an (n, m, p) matrix into a mesh, using marching_cubes.
 
@@ -103,6 +108,12 @@ def matrix_to_marching_cubes(matrix, pitch=1.0):
     -----------
     matrix : (n, m, p) bool
       Occupancy array
+    pitch : float or length-3 tuple of floats, optional
+      Voxel spacing in each dimension
+    threshold : float or None, optional
+      If specified, converts the input into a boolean
+      matrix by considering values above `threshold` as True
+
 
     Returns
     ----------
@@ -114,7 +125,10 @@ def matrix_to_marching_cubes(matrix, pitch=1.0):
 
     from ..base import Trimesh
 
-    matrix = np.asanyarray(matrix, dtype=bool)
+    if threshold is not None:
+        matrix = np.asanyarray(matrix) > threshold
+    else:
+        matrix = np.asanyarray(matrix, dtype=bool)
 
     rev_matrix = np.logical_not(matrix)  # Takes set about 0.
     # Add in padding so marching cubes can function properly with
diff --git a/trimesh/voxel/runlength.py b/trimesh/voxel/runlength.py
index 0b4a042c..f4c887ea 100644
--- a/trimesh/voxel/runlength.py
+++ b/trimesh/voxel/runlength.py
@@ -390,8 +390,9 @@ def sorted_rle_gather_1d(rle_data, ordered_indices):
                 start += next(data_iter)
             except StopIteration:
                 raise IndexError(
-                    "Index %d out of range of raw_values length %d" % (index, start)
+                    "Index %d out of range of raw_values length %d", index, start
                 )
+
         try:
             while index < start:
                 yield value
@@ -533,8 +534,9 @@ def sorted_brle_gather_1d(brle_data, ordered_indices):
                 start += next(data_iter)
             except StopIteration:
                 raise IndexError(
-                    "Index %d out of range of raw_values length %d" % (index, start)
+                    "Index %d out of range of raw_values length %d", index, start
                 )
+
         try:
             while index < start:
                 yield value


Test patch:
diff --git a/tests/test_iteration.py b/tests/test_iteration.py
new file mode 100644
index 00000000..917e5748
--- /dev/null
+++ b/tests/test_iteration.py
@@ -0,0 +1,65 @@
+from functools import reduce
+
+import numpy as np
+
+from trimesh.iteration import chain, reduce_cascade
+
+
+def test_reduce_cascade():
+    # the multiply will explode quickly past the integer maximum
+    def both(operation, items):
+        """
+        Run our cascaded reduce and regular reduce.
+        """
+
+        b = reduce_cascade(operation, items)
+
+        if len(items) > 0:
+            assert b == reduce(operation, items)
+
+        return b
+
+    for i in range(20):
+        data = np.arange(i)
+        c = both(items=data, operation=lambda a, b: a + b)
+
+        if i == 0:
+            assert c is None
+        else:
+            assert c == np.arange(i).sum()
+
+        # try a multiply
+        data = np.arange(i)
+        c = both(items=data, operation=lambda a, b: a * b)
+
+        if i == 0:
+            assert c is None
+        else:
+            assert c == np.prod(data)
+
+        # try a multiply
+        data = np.arange(i)[1:]
+        c = both(items=data, operation=lambda a, b: a * b)
+        if i <= 1:
+            assert c is None
+        else:
+            assert c == np.prod(data)
+
+    data = ["a", "b", "c", "d", "e", "f", "g"]
+    print("# reduce_pairwise\n-----------")
+    r = both(operation=lambda a, b: a + b, items=data)
+
+    assert r == "abcdefg"
+
+
+def test_chain():
+    # should work on iterables the same as `itertools.chain`
+    assert np.allclose(chain([1, 3], [4]), [1, 3, 4])
+    # should work with non-iterable single values
+    assert np.allclose(chain([1, 3], 4), [1, 3, 4])
+    # should filter out `None` arguments
+    assert np.allclose(chain([1, 3], None, 4, None), [1, 3, 4])
+
+
+if __name__ == "__main__":
+    test_reduce_cascade()
diff --git a/tests/test_remesh.py b/tests/test_remesh.py
index 2988a827..e5e3fc19 100644
--- a/tests/test_remesh.py
+++ b/tests/test_remesh.py
@@ -86,6 +86,10 @@ class SubDivideTest(g.unittest.TestCase):
         meshes = [g.trimesh.creation.box(), g.trimesh.creation.icosphere()]
 
         for m in meshes:
+            # set vertex positions as attributes for trivial check after subdivision
+            # make sure we're copying the array to avoid in-place check
+            m.vertex_attributes = {"pos": g.np.array(m.vertices) + 1.0}
+
             s = m.subdivide(face_index=[0, len(m.faces) - 1])
             # shouldn't have subdivided in-place
             assert len(s.faces) > len(m.faces)
@@ -93,6 +97,8 @@ class SubDivideTest(g.unittest.TestCase):
             assert g.np.isclose(m.area, s.area)
             # volume should be the same
             assert g.np.isclose(m.volume, s.volume)
+            # position attributes and actual vertices should be the same
+            assert g.np.allclose(s.vertex_attributes["pos"], s.vertices + 1.0)
 
             max_edge = m.scale / 50
             s = m.subdivide_to_size(max_edge=max_edge)
diff --git a/tests/test_util.py b/tests/test_util.py
index f8fbd2ef..8117bde7 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -75,16 +75,6 @@ class UtilTests(unittest.TestCase):
             # this is what should happen
             pass
 
-    def test_chain(self):
-        from trimesh.util import chain
-
-        # should work on iterables the same as `itertools.chain`
-        assert g.np.allclose(chain([1, 3], [4]), [1, 3, 4])
-        # should work with non-iterable single values
-        assert g.np.allclose(chain([1, 3], 4), [1, 3, 4])
-        # should filter out `None` arguments
-        assert g.np.allclose(chain([1, 3], None, 4, None), [1, 3, 4])
-
     def test_has_module(self):
         # built-in
         assert g.trimesh.util.has_module("collections")
diff --git a/tests/test_voxel.py b/tests/test_voxel.py
index 79a4bb52..f5d82c98 100644
--- a/tests/test_voxel.py
+++ b/tests/test_voxel.py
@@ -67,16 +67,29 @@ class VoxelGridTest(g.unittest.TestCase):
             g.log.warning("no skimage, skipping marching cubes test")
             return
 
+        march = g.trimesh.voxel.ops.matrix_to_marching_cubes
+
         # make sure offset is correct
         matrix = g.np.ones((3, 3, 3), dtype=bool)
-        mesh = g.trimesh.voxel.ops.matrix_to_marching_cubes(matrix=matrix)
+        mesh = march(matrix=matrix)
         assert mesh.is_watertight
 
-        mesh = g.trimesh.voxel.ops.matrix_to_marching_cubes(matrix=matrix).apply_scale(
-            3.0
-        )
+        mesh = march(matrix=matrix).apply_scale(3.0)
         assert mesh.is_watertight
 
+        # try an array full of a small number
+        matrix = g.np.full((3, 3, 3), 0.01, dtype=g.np.float64)
+        # set some to zero
+        matrix[:2, :2, :2] = 0.0
+
+        a = march(matrix)
+        assert a.is_watertight
+
+        # but above the threshold it should be not empty
+        b = march(matrix, threshold=-0.001)
+        assert b.is_watertight
+        assert b.volume > a.volume
+
     def test_marching_points(self):
         """
         Try marching cubes on points