Repository: mikedh/trimesh

Issue:
apply_transform and process have undocumented side effects
Hello,

Thanks for your library it's tremendously helpful.

I noticed that the `Trimesh.apply_transform` and `Trimesh.process` both have undocumented side effects.

The documentation for the `apply_transform` will be used here as an example of what I think is a problem. The function returns a `Trimesh` object : 

```
    def apply_transform(self, matrix):
        """
        Transform mesh by a homogeneous transformation matrix.
        Does the bookkeeping to avoid recomputing things so this function
        should be used rather than directly modifying self.vertices
        if possible.
        Parameters
        ------------
        matrix : (4, 4) float
          Homogeneous transformation matrix
        """
```

However, the apply_transform function not only returns a trimesh object, it also modifies the vertices of the called mesh, which is an undocumented side effect of the function.

As such, doing this : `mesh2 = mesh1.apply_transform(matrix)` also modified `mesh1`, which I think is a code smell and would yield many head scratching results (as it happened for me).

I would propose one of the following fixes in order to solve this problem : 

- Use an `inplace` keyword that would work as in the `pandas` library. When `inplace==True`, the caller is modified and no objects are returned, when `inplace==False` the caller isn't modified and a new instance is returned with the modification. I however understand that this would be a change that would break the standards of the library as it wouldn't be found anywhere else;
- `apply_transform` returns `None`. This way, it would be clear that `apply_transform` modifies the caller. This would be a breaking change, so I don't think this is a good solution. It would also break fluent interface;
- Finally, what I think would be the easiest solution is to add documentation to the function header to modify it to something like this : 

```
    def apply_transform(self, matrix):
        """
        Transform mesh by a homogeneous transformation matrix.
        Does the bookkeeping to avoid recomputing things so this function
        should be used rather than directly modifying self.vertices
        if possible.

        Parameters
        ------------
        matrix : (4, 4) float
          Homogeneous transformation matrix

        Returns
        -------
        Trimesh
          (self) A reference to the transformed instance mesh
        """
```

This last solution doesn't break any existing code, however it makes it very clear that `mesh1` is modified in the previous scenario. In order to figure out that `self` was returned, I needed to go in the source code. I believe this should be documented somewhere.

The same could be applied to the `process` function and to any other function that returns `self`.

Thanks

Discussion:

Patch:
diff --git a/Dockerfile b/Dockerfile
index e3f3548c..ceed0e0d 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -71,7 +71,6 @@ RUN pip install -e .[all]
 
 # check formatting
 RUN ruff trimesh
-RUN black --check trimesh
 
 # run pytest wrapped with xvfb for simple viewer tests
 RUN xvfb-run pytest --cov=trimesh \
diff --git a/docs/content/contributing.md b/docs/content/contributing.md
index 6eafd6d2..f8d48bde 100644
--- a/docs/content/contributing.md
+++ b/docs/content/contributing.md
@@ -71,9 +71,10 @@ if __name__ == '__main__':
 When you remove the embed and see the profile result you can then tweak the lines that are slow before finishing the function.
 
 ### Automatic Formatting
-The only check in that's required to pass in CI is `ruff`, which I usually run with:
+Trimesh uses `ruff` and `black` configured in `pyproject.toml`, you can run with:
 ```
 ruff . --fix
+black .
 ```
 It can fix a lot of formatting issues automatically. We also periodically run `black` to autoformat the codebase.
 
@@ -82,6 +83,7 @@ It can fix a lot of formatting issues automatically. We also periodically run `b
 
 Trimesh uses the [Sphinx Numpy-style](https://www.sphinx-doc.org/en/master/usage/extensions/example_numpy.html#example-numpy) docstrings which get parsed into the API reference page. 
 
+
 ## General Tips
 
 Python can be fast but only when you use it as little as possible. In general, if you ever have a block which loops through faces and vertices it will be basically unusable with even moderately sized meshes. All operations on face or vertex arrays should be vectorized numpy operations unless absolutely unavoidable. Profiling helps figure out what is slow, but some general advice:
diff --git a/docs/content/install.md b/docs/content/install.md
index 9a769722..18c3a3cc 100644
--- a/docs/content/install.md
+++ b/docs/content/install.md
@@ -24,8 +24,7 @@ pip install trimesh[all]
 ```
 
 
-Conda Packages
---------------
+## Conda Packages
 
 If you prefer a `conda` environment, `trimesh` is available on `conda-forge` ([trimesh-feedstock repo](https://github.com/conda-forge/trimesh-feedstock))
 
@@ -34,20 +33,11 @@ If you install [Miniconda](https://docs.conda.io/projects/miniconda/en/latest/)
 ```
 conda install -c conda-forge trimesh
 ```
-      
-Ubuntu-Debian Notes
--------------------
 
-Blender and openSCAD are soft dependencies used for boolean operations with subprocess, you can get them with `apt`:
-
-```
-sudo apt-get install blender
-```
-
-Dependency Overview
+## Dependency Overview
 --------------------
 
-Trimesh has a lot of soft-required upstream packages. We try to make sure they're active and big-ish. Here's a quick summary of what they're used for.
+Trimesh has a lot of soft-required upstream packages, and we try to make sure they're actively maintained. Here's a quick summary of what they're used for:
 
  
 | Package | Description | Alternatives | Level |
@@ -81,3 +71,16 @@ Trimesh has a lot of soft-required upstream packages. We try to make sure they'r
 |`pytest-cov`| A plugin to calculate test coverage. | | `test`|
 |`pyinstrument`| A sampling based profiler for performance tweaking. | | `test`|
 |`vhacdx`| A binding for VHACD which provides convex decompositions | | `recommend`|
+
+## Adding A Dependency
+
+If there's no way to implement something reasonably in vectorized Python or there is a mature minimal C++ or Rust implementation of something useful and complicated we may add a dependency. If it's a major, active project with few dependencies (i.e. `jinja2`) that's probably fine. Otherwise it's a lot more of a commitment than just implementing the function in Python however. An example of this is `embree`, Intel's ray check engine: it is a super complicated thing to do well and 50-100x faster than Python ray checks.
+
+There are a few projects that we've forked into the [`trimesh`](https://github.com/trimesh/) GitHub organization which you can take a look at. The general idea of the requirements for a new compiled dependency are:
+
+- is actively maintained and has an MIT/BSD compatible license.
+- has all source code in the repository or as a submodule, i.e. no mysterious binary blobs.
+- binding preferably uses [pybind11](https://pybind11.readthedocs.io/en/stable/index.html), [nanobind](https://github.com/wjakob/nanobind) or [maturin/py03](https://github.com/PyO3/maturin) for Rust projects. Cython is also OK but other options are preferable if possible. 
+- uses `cibuildwheel` to publish releases configured in `pyproject.toml`.
+- has unit tests which run in CI
+- has minimal dependencies: ideally only `numpy`.
\ No newline at end of file
diff --git a/examples/offscreen_render.py b/examples/offscreen_render.py
index 797ce8b9..a8db740a 100644
--- a/examples/offscreen_render.py
+++ b/examples/offscreen_render.py
@@ -17,11 +17,11 @@ if __name__ == "__main__":
     # a 45 degree homogeneous rotation matrix around
     # the Y axis at the scene centroid
     rotate = trimesh.transformations.rotation_matrix(
-        angle=np.radians(10.0), direction=[0, 1, 0], point=scene.centroid
+        angle=np.radians(30.0), direction=[1, 0, 0], point=scene.centroid
     )
 
-    for i in range(4):
-        trimesh.constants.log.info("Saving image %d", i)
+    for i in range(10):
+        trimesh.constants.log.info(f"Saving image {i}")
 
         # rotate the camera view transform
         camera_old, _geometry = scene.graph[scene.camera.name]
diff --git a/pyproject.toml b/pyproject.toml
index 07a1979c..2b4aa8c3 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,7 +5,7 @@ requires = ["setuptools >= 61.0", "wheel"]
 [project]
 name = "trimesh"
 requires-python = ">=3.7"
-version = "4.0.0"
+version = "4.0.1"
 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/base.py b/trimesh/base.py
index 8ef060d1..f8f9ca2e 100644
--- a/trimesh/base.py
+++ b/trimesh/base.py
@@ -37,7 +37,7 @@ from . import (
     units,
     util,
 )
-from .constants import log, log_time, tol
+from .constants import log, tol
 from .exceptions import ExceptionWrapper
 from .exchange.export import export_mesh
 from .parent import Geometry3D
@@ -277,7 +277,7 @@ class Trimesh(Geometry3D):
         return self._data.mutable
 
     @mutable.setter
-    def mutable(self, value: bool):
+    def mutable(self, value: bool) -> None:
         """
         Set the mutability of the current mesh.
 
@@ -342,7 +342,7 @@ class Trimesh(Geometry3D):
         return sparse
 
     @property
-    def face_normals(self):
+    def face_normals(self) -> NDArray[float64]:
         """
         Return the unit normal vector for each face.
 
@@ -394,7 +394,7 @@ class Trimesh(Geometry3D):
         return padded
 
     @face_normals.setter
-    def face_normals(self, values):
+    def face_normals(self, values: NDArray[float64]) -> None:
         """
         Assign values to face normals.
 
@@ -434,7 +434,7 @@ class Trimesh(Geometry3D):
         self._cache["face_normals"] = values
 
     @property
-    def vertices(self):
+    def vertices(self) -> NDArray[float64]:
         """
         The vertices of the mesh.
 
@@ -451,7 +451,7 @@ class Trimesh(Geometry3D):
         return self._data.get("vertices", np.empty(shape=(0, 3), dtype=float64))
 
     @vertices.setter
-    def vertices(self, values):
+    def vertices(self, values: NDArray[float64]):
         """
         Assign vertex values to the mesh.
 
@@ -465,7 +465,7 @@ class Trimesh(Geometry3D):
         self._data["vertices"] = np.asanyarray(values, order="C", dtype=float64)
 
     @caching.cache_decorator
-    def vertex_normals(self):
+    def vertex_normals(self) -> NDArray[float64]:
         """
         The vertex normals of the mesh. If the normals were loaded
         we check to make sure we have the same number of vertex
@@ -491,7 +491,7 @@ class Trimesh(Geometry3D):
         return vertex_normals
 
     @vertex_normals.setter
-    def vertex_normals(self, values: NDArray[float64]):
+    def vertex_normals(self, values: NDArray[float64]) -> None:
         """
         Assign values to vertex normals.
 
@@ -608,7 +608,7 @@ class Trimesh(Geometry3D):
         return centroid
 
     @property
-    def center_mass(self):
+    def center_mass(self) -> NDArray[float64]:
         """
         The point in space which is the center of mass/volume.
 
@@ -620,7 +620,7 @@ class Trimesh(Geometry3D):
         return self.mass_properties.center_mass
 
     @center_mass.setter
-    def center_mass(self, value):
+    def center_mass(self, value: NDArray[float64]) -> None:
         """
         Override the point in space which is the center of mass and volume.
 
@@ -648,7 +648,7 @@ class Trimesh(Geometry3D):
         return self.mass_properties.density
 
     @density.setter
-    def density(self, value: float):
+    def density(self, value: float) -> None:
         """
         Set the density of the primitive.
 
@@ -1138,9 +1138,9 @@ class Trimesh(Geometry3D):
         self,
         merge_tex: Optional[bool] = None,
         merge_norm: Optional[bool] = None,
-        digits_vertex: None = None,
-        digits_norm: None = None,
-        digits_uv: None = None,
+        digits_vertex: Optional[bool] = None,
+        digits_norm: Optional[bool] = None,
+        digits_uv: Optional[bool] = None,
     ) -> None:
         """
         Removes duplicate vertices grouped by position and
@@ -1336,7 +1336,7 @@ class Trimesh(Geometry3D):
         )
         self.update_faces(self.unique_faces())
 
-    def rezero(self):
+    def rezero(self) -> None:
         """
         Translate the mesh so that all vertex vertices are positive.
 
@@ -1709,7 +1709,7 @@ class Trimesh(Geometry3D):
         )
         self.update_faces(self.nondegenerate_faces(height=height))
 
-    def nondegenerate_faces(self, height=tol.merge) -> NDArray[bool]:
+    def nondegenerate_faces(self, height: float = tol.merge) -> NDArray[bool]:
         """
         Remove degenerate faces (faces without 3 unique vertex indices)
         from the current mesh.
@@ -1855,7 +1855,7 @@ class Trimesh(Geometry3D):
 
         return on_hull
 
-    def fix_normals(self, multibody: Optional[bool] = None):
+    def fix_normals(self, multibody: Optional[bool] = None) -> None:
         """
         Find and fix problems with self.face_normals and self.faces
         winding direction.
@@ -1885,7 +1885,7 @@ class Trimesh(Geometry3D):
         """
         return repair.fill_holes(self)
 
-    def register(self, other, **kwargs):
+    def register(self, other: Geometry3D, **kwargs):
         """
         Align a mesh with another mesh or a PointCloud using
         the principal axes of inertia as a starting point which
@@ -1917,7 +1917,11 @@ class Trimesh(Geometry3D):
         return mesh_to_other, cost
 
     def compute_stable_poses(
-        self, center_mass=None, sigma=0.0, n_samples=1, threshold=0.0
+        self,
+        center_mass: Optional[NDArray[float64]] = None,
+        sigma: float = 0.0,
+        n_samples: int = 1,
+        threshold: float = 0.0,
     ):
         """
         Computes stable orientations of a mesh and their quasi-static probabilities.
@@ -2112,38 +2116,48 @@ class Trimesh(Geometry3D):
         result = Trimesh(vertices=new_vertices, faces=new_faces, process=False)
         return result
 
-    @log_time
     def smoothed(self, **kwargs):
         """
-        Return a version of the current mesh which will render
-        nicely, without changing source mesh.
+        DEPRECATED: use `mesh.smooth_shaded` or `trimesh.graph.smooth_shade(mesh)`
+        """
+        warnings.warn(
+            "`mesh.smoothed()` is deprected and will be removed in March 2024: "
+            + "use `mesh.smooth_shaded` or `trimesh.graph.smooth_shade(mesh)`",
+            category=DeprecationWarning,
+            stacklevel=2,
+        )
+        # run smoothing
+        return self.smooth_shaded
 
-        Parameters
-        -------------
-        angle : float or None
-          Angle in radians face pairs with angles
-          smaller than this will appear smoothed
-        facet_minarea : float or None
-          Minimum area fraction to consider
-          IE for `facets_minarea=25` only facets larger
-          than `mesh.area / 25` will be considered.
+    @property
+    def smooth_shaded(self):
+        """
+        Smooth shading in OpenGL relies on which vertices are shared,
+        this function will disconnect regions above an angle threshold
+        and return a non-watertight version which will look better
+        in an OpenGL rendering context.
+
+        If you would like to use non-default arguments see `graph.smooth_shade`.
 
         Returns
         ---------
-        smoothed : trimesh.Trimesh
-          Non watertight version of current mesh
-          which will render nicely with smooth shading
+        smooth_shaded : trimesh.Trimesh
+          Non watertight version of current mesh.
         """
-
-        # smooth should be recomputed if visuals change
+        # key this also by the visual properties
+        # but store it in the mesh cache
         self.visual._verify_hash()
-        cached = self.visual._cache["smoothed"]
-        if cached is not None:
-            return cached
-        # run smoothing
-        smoothed = graph.smoothed(self, **kwargs)
-        self.visual._cache["smoothed"] = smoothed
-        return smoothed
+
+        cache = self.visual._cache
+        # needs to be dumped whenever visual or mesh changes
+        key = f"smooth_shaded_{hash(self.visual)}_{hash(self)}"
+        if key in cache:
+            return cache[key]
+        smooth = graph.smooth_shade(self)
+
+        # store it in the mesh cache which dumps when vertices change
+        cache[key] = smooth
+        return smooth
 
     @property
     def visual(self):
@@ -2216,7 +2230,12 @@ class Trimesh(Geometry3D):
 
         return path
 
-    def section_multiplane(self, plane_origin, plane_normal, heights):
+    def section_multiplane(
+        self,
+        plane_origin: NDArray[float64],
+        plane_normal: NDArray[float64],
+        heights: NDArray[float64],
+    ):
         """
         Return multiple parallel cross sections of the current
         mesh in 2D.
diff --git a/trimesh/graph.py b/trimesh/graph.py
index a816402d..a9160d73 100644
--- a/trimesh/graph.py
+++ b/trimesh/graph.py
@@ -9,12 +9,14 @@ Currently uses networkx or scipy.sparse.csgraph backend.
 """
 
 import collections
+import warnings
 
 import numpy as np
 
 from . import exceptions, grouping, util
 from .constants import log, tol
 from .geometry import faces_to_edges
+from .typed import Optional
 
 try:
     from scipy.sparse import coo_matrix, csgraph
@@ -741,9 +743,24 @@ def neighbors(edges, max_index=None, directed=False):
     return array
 
 
-def smoothed(mesh, angle=None, facet_minarea=10):
+def smoothed(*args, **kwargs):
     """
-    Return a non- watertight version of the mesh which
+    DEPRECATED: use `trimesh.graph.smooth_shade(mesh, ...)`
+    """
+    warnings.warn(
+        "`trimesh.graph.smoothed` is deprected and will be removed in March 2024: "
+        + "use `trimesh.graph.smooth_shade(mesh, ...)`",
+        category=DeprecationWarning,
+        stacklevel=2,
+    )
+    return smooth_shade(*args, **kwargs)
+
+
+def smooth_shade(
+    mesh, angle: Optional[float] = None, facet_minarea: Optional[float] = 10.0
+):
+    """
+    Return a non-watertight version of the mesh which
     will render nicely with smooth shading by
     disconnecting faces at sharp angles to each other.
 
diff --git a/trimesh/grouping.py b/trimesh/grouping.py
index 44a6c12c..36abcfca 100644
--- a/trimesh/grouping.py
+++ b/trimesh/grouping.py
@@ -781,6 +781,10 @@ def blocks(data, min_len=2, max_len=np.inf, wrap=False, digits=None, only_nonzer
         if only_nonzero and not bool(data[0]):
             return blocks
 
+        # if all values are True or False we can exit
+        if len(blocks) == 1 and len(blocks[0]) == len(data):
+            return blocks
+
         # so now first point equals last point, so the cases are:
         # - first and last point are in a block: combine two blocks
         # - first OR last point are in block: add other point to block
diff --git a/trimesh/interval.py b/trimesh/interval.py
index 34fa856e..73782d5b 100644
--- a/trimesh/interval.py
+++ b/trimesh/interval.py
@@ -8,114 +8,100 @@ Deal with 1D intervals which are defined by:
 
 import numpy as np
 
+from .typed import NDArray, float64
 
-def check(a, b, digits):
+
+def intersection(a: NDArray[float64], b: NDArray[float64]) -> NDArray[float64]:
     """
-    Check input ranges, convert them to vector form,
-    and get a fixed precision integer version of them.
+    Given pairs of ranges merge them in to
+    one range if they overlap.
 
     Parameters
     --------------
-    a : (2, ) or (2, n) float
+    a : (2, ) or (n, 2)
       Start and end of a 1D interval
-    b : (2, ) or (2, n) float
+    b : (2, ) float
       Start and end of a 1D interval
-    digits : int
-      How many digits to consider
 
     Returns
     --------------
-    a : (2, n) float
-      Ranges as vector
-    b : (2, n) float
-      Ranges as vector
-    a_int : (2, n) int64
-      Ranges rounded to digits, as vector
-    b_int : (2, n) int64
-      Ranges rounded to digits, as vector
-    is_1D : bool
-      If True, input was single pair of ranges
+    inter : (2, ) or (2, 2) float
+      The unioned range from the two inputs,
+      if not `inter.ptp(axis=1)` will be zero.
     """
     a = np.array(a, dtype=np.float64)
     b = np.array(b, dtype=np.float64)
 
-    if a.shape != b.shape or a.shape[-1] != 2:
-        raise ValueError("ranges must be identical and (2,)!")
-
-    # if input was single interval reshape it here
-    is_1D = False
-    if len(a.shape) == 1:
-        a = a.reshape((-1, 2))
-        b = b.reshape((-1, 2))
-        is_1D = True
+    # convert to vectorized form
+    is_1D = a.shape == (2,)
+    a = a.reshape((-1, 2))
+    b = b.reshape((-1, 2))
 
-    # make sure ranges are sorted
+    # make sure they're min-max
     a.sort(axis=1)
     b.sort(axis=1)
+    a_low, a_high = a.T
+    b_low, b_high = b.T
 
-    # compare in fixed point as integers
-    a_int = (a * 10**digits).round().astype(np.int64)
-    b_int = (b * 10**digits).round().astype(np.int64)
-
-    return a, b, a_int, b_int, is_1D
+    # do the checks
+    check = np.logical_not(np.logical_or(b_low >= a_high, a_low >= b_high))
+    overlap = np.zeros(a.shape, dtype=np.float64)
+    overlap[check] = np.column_stack(
+        (
+            np.array([a_low[check], b_low[check]]).max(axis=0),
+            np.array([a_high[check], b_high[check]]).min(axis=0),
+        )
+    )
 
+    if is_1D:
+        return overlap[0]
 
-def intersection(a, b, digits=8):
-    """
-    Given a pair of ranges, merge them in to
-    one range if they overlap at all
+    return overlap
 
-    Parameters
-    --------------
-    a : (2, ) float
-      Start and end of a 1D interval
-    b : (2, ) float
-      Start and end of a 1D interval
-    digits : int
-      How many digits to consider
 
-    Returns
-    --------------
-    intersects : bool or (n,) bool
-      Indicates if the ranges overlap at all
-    new_range : (2, ) or (2, 2) float
-      The unioned range from the two inputs,
-      or both of the original ranges if not overlapping
+def union(intervals: NDArray[float64], sort: bool = True) -> NDArray[float64]:
     """
-    # check shape and convert
-    a, b, a_int, b_int, is_1D = check(a, b, digits)
-
-    # what are the starting and ending points of the overlap
-    overlap = np.zeros(a.shape, dtype=np.float64)
-
-    # A fully overlaps B
-    current = np.logical_and(a_int[:, 0] <= b_int[:, 0], a_int[:, 1] >= b_int[:, 1])
-    overlap[current] = b[current]
-
-    # B fully overlaps A
-    current = np.logical_and(a_int[:, 0] >= b_int[:, 0], a_int[:, 1] <= b_int[:, 1])
-    overlap[current] = a[current]
+    For array of multiple intervals union them all into
+    the subset of intervals.
 
-    # A starts B ends
-    # A:, 0   B:, 0     A:, 1        B:, 1
-    current = np.logical_and(
-        np.logical_and(a_int[:, 0] <= b_int[:, 0], b_int[:, 0] < a_int[:, 1]),
-        a_int[:, 1] < b_int[:, 1],
-    )
-    overlap[current] = np.column_stack([b[current][:, 0], a[current][:, 1]])
-
-    # B starts A ends
-    # B:, 0  A:, 0    B:, 1  A:, 1
-    current = np.logical_and(
-        np.logical_and(b_int[:, 0] <= a_int[:, 0], a_int[:, 0] < b_int[:, 1]),
-        b_int[:, 1] < a_int[:, 1],
-    )
-    overlap[current] = np.column_stack([a[current][:, 0], b[current][:, 1]])
+    For example:
+    `intervals = [[1,2], [2,3]] -> [[1, 3]]`
+    `intervals = [[1,2], [2.5,3]] -> [[1, 2], [2.5, 3]]`
 
-    # is range overlapping at all
-    intersects = overlap.ptp(axis=1) > 10**-digits
 
-    if is_1D:
-        return intersects[0], overlap[0]
+    Parameters
+    ------------
+    intervals : (n, 2)
+      Pairs of `(min, max)` values.
+    sort
+      If the array is already ordered into (min, max) pairs
+      and then pairs sorted by minimum value you can skip the
+      sorting in this function.
 
-    return intersects, overlap
+    Returns
+    ----------
+    unioned : (m, 2)
+      New intervals where `m <= n`
+    """
+    if len(intervals) == 0:
+        return np.empty(0)
+
+    # if the intervals have not been pre-sorted we should apply our sorting logic
+    # you would only skip this if you are subsetting a larger list elsewhere.
+    if sort:
+        # copy inputs and make sure they are (min, max) pairs
+        intervals = np.sort(intervals, axis=1)
+        # order them by lowest starting point
+        intervals = intervals[intervals[:, 0].argsort()]
+
+    # we know we will have at least one interval
+    unions = [intervals[0]]
+
+    for begin, end in intervals[1:]:
+        if unions[-1][1] >= begin:
+            #
+            unions[-1][1] = max(unions[-1][1], end)
+        else:
+            unions.append([begin, end])
+
+    return np.array(unions)
diff --git a/trimesh/path/segments.py b/trimesh/path/segments.py
index 7ab3ebde..566be34a 100644
--- a/trimesh/path/segments.py
+++ b/trimesh/path/segments.py
@@ -7,11 +7,13 @@ Deal with (n, 2, 3) line segments.
 
 import numpy as np
 
-from .. import geometry, grouping, interval, transformations, util
+from .. import geometry, interval, transformations, util
 from ..constants import tol
+from ..grouping import group_rows, unique_rows
+from ..typed import NDArray, float64
 
 
-def segments_to_parameters(segments):
+def segments_to_parameters(segments: NDArray[float64]):
     """
     For 3D line segments defined by two points, turn
     them in to an origin defined as the closest point along
@@ -56,7 +58,9 @@ def segments_to_parameters(segments):
     return origins, vectors, parameters
 
 
-def parameters_to_segments(origins, vectors, parameters):
+def parameters_to_segments(
+    origins: NDArray[float64], vectors: NDArray[float64], parameters: NDArray[float64]
+):
     """
     Convert a parametric line segment representation to
     a two point line segment representation
@@ -146,6 +150,44 @@ def colinear_pairs(segments, radius=0.01, angle=0.01, length=None):
     return colinear
 
 
+def clean(segments: NDArray[float64], digits: int = 10) -> NDArray[float64]:
+    """
+    Clean up line segments by unioning the ranges of colinear segments.
+
+    Parameters
+    ------------
+    segments : (n, 2, 2) or (n, 2, 3)
+      Line segments in space.
+    digits
+      How many digits to consider.
+
+    Returns
+    -----------
+    cleaned : (m, 2, 2) or (m, 2, 3)
+      Where `m <= n`
+    """
+    # convert segments to parameterized origins
+    # which are the closest point on the line to
+    # the actual zero- origin
+    origins, vectors, param = segments_to_parameters(segments)
+
+    # make sure parameters are in min-max order
+    param.sort(axis=1)
+
+    # collect new unified paramameters
+    p, o, v = [], [], []
+    for g in group_rows(np.column_stack((origins, vectors)), digits=digits):
+        # union the intervals sorting ourselves to skip the `sort(axis=1)` we did above
+        group = param[g]
+        u = interval.union(group[group[:, 0].argsort()], sort=False)
+        p.extend(u)
+        # use the origins for the subsetted union
+        o.extend(origins[g[: len(u)]])
+        v.extend(vectors[g[: len(u)]])
+
+    return parameters_to_segments(o, v, p)
+
+
 def split(segments, points, atol=1e-5):
     """
     Find any points that lie on a segment (not an endpoint)
@@ -231,15 +273,15 @@ def unique(segments, digits=5):
     segments = np.asanyarray(segments, dtype=np.float64)
 
     # find segments as unique indexes so we can find duplicates
-    inverse = grouping.unique_rows(
-        segments.reshape((-1, segments.shape[2])), digits=digits
-    )[1].reshape((-1, 2))
+    inverse = unique_rows(segments.reshape((-1, segments.shape[2])), digits=digits)[
+        1
+    ].reshape((-1, 2))
     # make sure rows are sorted
     inverse.sort(axis=1)
     # remove segments where both indexes are the same
     mask = np.zeros(len(segments), dtype=bool)
     # only include the first occurrence of a segment
-    mask[grouping.unique_rows(inverse)[0]] = True
+    mask[unique_rows(inverse)[0]] = True
     # remove segments that are zero-length
     mask[inverse[:, 0] == inverse[:, 1]] = False
     # apply the unique mask
@@ -248,68 +290,6 @@ def unique(segments, digits=5):
     return unique
 
 
-def overlap(origins, vectors, params):
-    """
-    Find the overlap of two parallel line segments.
-
-    Parameters
-    ------------
-    origins : (2, 3) float
-       Origin points of lines in space
-    vectors : (2, 3) float
-       Unit direction vectors of lines
-    params : (2, 2) float
-       Two (start, end) distance pairs
-
-    Returns
-    ------------
-    length : float
-       Overlapping length
-    overlap : (n, 2, 3) float
-       Line segments for overlapping distance
-    """
-    # copy inputs and make sure shape is correct
-    origins = np.array(origins).reshape((2, 3))
-    vectors = np.array(vectors).reshape((2, 3))
-    params = np.array(params).reshape((2, 2))
-
-    if tol.strict:
-        # convert input to parameters before flipping
-        # to make sure we didn't screw it up
-        truth = parameters_to_segments(origins, vectors, params)
-
-    # this function only works on parallel lines
-    dot = np.dot(*vectors)
-    assert np.isclose(np.abs(dot), 1.0, atol=0.01)
-
-    # if two vectors are reversed
-    if dot < 0.0:
-        # reverse direction vector
-        vectors[1] *= -1.0
-        # negate parameters
-        params[1] *= -1.0
-
-    if tol.strict:
-        # do a check to make sure our reversal didn't
-        # inadvertently give us incorrect segments
-        assert np.allclose(truth, parameters_to_segments(origins, vectors, params))
-
-    # merge the parameter ranges
-    ok, new_range = interval.intersection(*params)
-
-    if not ok:
-        return 0.0, np.array([])
-
-    # create the overlapping segment pairs (2, 2, 3)
-    segments = np.array(
-        [o + v * new_range.reshape((-1, 1)) for o, v in zip(origins, vectors)]
-    )
-    # get the length of the new range
-    length = new_range.ptp()
-
-    return length, segments
-
-
 def extrude(segments, height, double_sided=False):
     """
     Extrude 2D line segments into 3D triangles.
@@ -336,11 +316,16 @@ def extrude(segments, height, double_sided=False):
 
     # we are creating two vertices  triangles for every 2D line segment
     # on the segments of the 2D triangulation
-    vertices = np.tile(segments.reshape((-1, 2)), 2).reshape((-1, 2))
-    vertices = np.column_stack((vertices, np.tile([0, height, 0, height], len(segments))))
-    faces = np.tile([3, 1, 2, 2, 1, 0], (len(segments), 1))
-    faces += np.arange(len(segments)).reshape((-1, 1)) * 4
-    faces = faces.reshape((-1, 3))
+    vertices = np.column_stack(
+        (
+            np.tile(segments.reshape((-1, 2)), 2).reshape((-1, 2)),
+            np.tile([0, height, 0, height], len(segments)),
+        )
+    )
+    faces = (
+        np.tile([3, 1, 2, 2, 1, 0], (len(segments), 1))
+        + np.arange(len(segments)).reshape((-1, 1)) * 4
+    ).reshape((-1, 3))
 
     if double_sided:
         # stack so they will render from the back
@@ -383,7 +368,7 @@ def resample(segments, maxlen, return_index=False, return_count=False):
 
     Parameters
     -------------
-    segments : (n, 2, 2) float
+    segments : (n, 2, 2|3) float
       2D line segments
     maxlen : float
       The maximum length of a line segment
@@ -394,7 +379,7 @@ def resample(segments, maxlen, return_index=False, return_count=False):
 
     Returns
     -------------
-    resampled : (m, 2, 3) float
+    resampled : (m, 2, 2|3) float
       Line segments where no segment is longer than maxlen
     index : (m,) int
       [OPTIONAL] The index of segments resampled came from
@@ -404,6 +389,10 @@ def resample(segments, maxlen, return_index=False, return_count=False):
     # check arguments
     maxlen = float(maxlen)
     segments = np.array(segments, dtype=np.float64)
+    if len(segments.shape) != 3:
+        raise ValueError(f"{segments.shape} != (n, 2, 2|3)")
+
+    dimension = segments.shape[2]
 
     # shortcut for endpoints
     pt1 = segments[:, 0]
@@ -429,13 +418,13 @@ def resample(segments, maxlen, return_index=False, return_count=False):
         # the vector for each incremental length
         increment = vec[mask] / split
         # stack the increment vector into the shape needed
-        v = tile(increment, split + 1).reshape((-1, 3)) * tile(
+        v = tile(increment, split + 1).reshape((-1, dimension)) * tile(
             np.arange(split + 1), len(increment)
         ).reshape((-1, 1))
         # stack the origin points correctly
-        o = tile(pt1[mask], split + 1).reshape((-1, 3))
+        o = tile(pt1[mask], split + 1).reshape((-1, dimension))
         # now get each segment as an (split, 3) polyline
-        poly = (o + v).reshape((-1, split + 1, 3))
+        poly = (o + v).reshape((-1, split + 1, dimension))
         # save the resulting segments
         # magical slicing is equivalent to:
         # > [p[stack] for p in poly]
diff --git a/trimesh/rendering.py b/trimesh/rendering.py
index 8aea52f9..1f81f882 100644
--- a/trimesh/rendering.py
+++ b/trimesh/rendering.py
@@ -69,6 +69,7 @@ def mesh_to_vertexlist(mesh, group=None, smooth=True, smooth_threshold=60000):
     --------------
     args : (7,) tuple
       Args for vertex list constructor
+
     """
 
     if hasattr(mesh.visual, "uv"):
@@ -107,7 +108,7 @@ def mesh_to_vertexlist(mesh, group=None, smooth=True, smooth_threshold=60000):
         # if we have a small number of faces and colors defined
         # smooth the  mesh by merging vertices of faces below
         # the threshold angle
-        mesh = mesh.smoothed()
+        mesh = mesh.smooth_shaded
         vertex_count = len(mesh.vertices)
         normals = mesh.vertex_normals.reshape(-1).tolist()
         faces = mesh.faces.reshape(-1).tolist()
@@ -134,6 +135,7 @@ def mesh_to_vertexlist(mesh, group=None, smooth=True, smooth_threshold=60000):
         ("n3f/static", normals),
         color_gl,
     )
+
     return args
 
 
diff --git a/trimesh/scene/scene.py b/trimesh/scene/scene.py
index d9326763..1fb80d99 100644
--- a/trimesh/scene/scene.py
+++ b/trimesh/scene/scene.py
@@ -252,6 +252,7 @@ class Scene(Geometry3D):
         """
         # avoid accessing attribute in tight loop
         geometry = self.geometry
+        # hash of geometry and transforms
         # start with the last modified time of the scene graph
         hashable = [hex(self.graph.transforms.__hash__())]
         # take the re-hex string of the hash
diff --git a/trimesh/viewer/windowed.py b/trimesh/viewer/windowed.py
index 7a360f8f..32704450 100644
--- a/trimesh/viewer/windowed.py
+++ b/trimesh/viewer/windowed.py
@@ -229,7 +229,7 @@ class SceneViewer(pyglet.window.Window):
         for name, geom in self.scene.geometry.items():
             if geom.is_empty:
                 continue
-            if geometry_hash(geom) == self.vertex_list_hash.get(name):
+            if _geometry_hash(geom) == self.vertex_list_hash.get(name):
                 continue
             self.add_geometry(name=name, geometry=geom, smooth=bool(self._smooth))
 
@@ -263,7 +263,7 @@ class SceneViewer(pyglet.window.Window):
         # create the indexed vertex list
         self.vertex_list[name] = self.batch.add_indexed(*args)
         # save the hash of the geometry
-        self.vertex_list_hash[name] = geometry_hash(geometry)
+        self.vertex_list_hash[name] = _geometry_hash(geometry)
         # save the rendering mode from the constructor args
         self.vertex_list_mode[name] = args[1]
 
@@ -731,7 +731,6 @@ class SceneViewer(pyglet.window.Window):
 
             # get the transform from world to geometry and mesh name
             transform, geometry_name = graph.get(current_node)
-
             # if no geometry at this frame continue without rendering
             if geometry_name is None or geometry_name not in self.vertex_list_mode:
                 continue
@@ -749,7 +748,6 @@ class SceneViewer(pyglet.window.Window):
             mesh = geometry[geometry_name]
             if mesh.is_empty:
                 continue
-
             # get the GL mode of the current geometry
             mode = self.vertex_list_mode[geometry_name]
 
@@ -842,7 +840,7 @@ class SceneViewer(pyglet.window.Window):
         return file_obj
 
 
-def geometry_hash(geometry):
+def _geometry_hash(geometry):
     """
     Get a hash for a geometry object
 


Test patch:
diff --git a/tests/regression.py b/tests/regression.py
index e554dec9..587a2436 100644
--- a/tests/regression.py
+++ b/tests/regression.py
@@ -31,7 +31,7 @@ def typical_application():
 
         faces = mesh.facets[mesh.facets_area.argmax()]
         outline = mesh.outline(faces)  # NOQA
-        smoothed = mesh.smoothed()  # NOQA
+        smoothed = mesh.smooth_shaded  # NOQA
 
         assert mesh.volume > 0.0
 
diff --git a/tests/test_cache.py b/tests/test_cache.py
index 6e515c5b..d4bc4715 100644
--- a/tests/test_cache.py
+++ b/tests/test_cache.py
@@ -356,6 +356,24 @@ class CacheTest(g.unittest.TestCase):
         m.process(validate=True)
         assert m.triangles.shape == (1, 3, 3)
 
+    def test_smooth_shade(self, count=10):
+        # test to make sure the smooth shaded copy is cached correctly
+        mesh = g.trimesh.creation.cylinder(radius=1, height=10)
+        scene = g.trimesh.Scene({"mesh": mesh})
+
+        initial = scene.camera_transform.copy()
+
+        hashes = set()
+        for n in range(count):
+            angle = g.np.pi * n / count
+            matrix = g.trimesh.transformations.rotation_matrix(angle, [1, 0, 0])
+            scene.geometry["mesh"].apply_transform(matrix)
+            hashes.add(hash(scene.geometry["mesh"].smooth_shaded))
+            scene.camera_transform = initial
+
+        # the smooth shade should be unique for every transform
+        assert len(hashes) == count
+
 
 if __name__ == "__main__":
     g.trimesh.util.attach_to_log()
diff --git a/tests/test_color.py b/tests/test_color.py
index 88f29387..83e95baa 100644
--- a/tests/test_color.py
+++ b/tests/test_color.py
@@ -123,7 +123,7 @@ class VisualTest(g.unittest.TestCase):
         m = g.get_mesh("featuretype.STL")
 
         # will put smoothed mesh into visuals cache
-        s = m.smoothed()
+        s = m.smooth_shaded
         # every color should be default color
         assert s.visual.face_colors.ptp(axis=0).max() == 0
 
@@ -131,16 +131,16 @@ class VisualTest(g.unittest.TestCase):
         m.visual.face_colors[0] = [255, 0, 0, 255]
 
         # cache should be dumped yo
-        s1 = m.smoothed()
+        s1 = m.smooth_shaded
         assert s1.visual.face_colors.ptp(axis=0).max() != 0
 
         # do the same check on vertex color
         m = g.get_mesh("featuretype.STL")
-        s = m.smoothed()
+        s = m.smooth_shaded
         # every color should be default color
         assert s.visual.vertex_colors.ptp(axis=0).max() == 0
         m.visual.vertex_colors[g.np.arange(10)] = [255, 0, 0, 255]
-        s1 = m.smoothed()
+        s1 = m.smooth_shaded
         assert s1.visual.face_colors.ptp(axis=0).max() != 0
 
     def test_vertex(self):
diff --git a/tests/test_graph.py b/tests/test_graph.py
index a389e238..035e3f10 100644
--- a/tests/test_graph.py
+++ b/tests/test_graph.py
@@ -113,7 +113,7 @@ class GraphTest(g.unittest.TestCase):
 
         for name in ["ADIS16480.STL", "featuretype.STL"]:
             mesh = g.get_mesh(name)
-            assert len(mesh.faces) == len(mesh.smoothed().faces)
+            assert len(mesh.faces) == len(mesh.smooth_shaded.faces)
 
     def test_engines(self):
         edges = g.np.arange(10).reshape((-1, 2))
diff --git a/tests/test_grouping.py b/tests/test_grouping.py
index 583d4bde..52a75d16 100644
--- a/tests/test_grouping.py
+++ b/tests/test_grouping.py
@@ -78,6 +78,22 @@ class GroupTests(g.unittest.TestCase):
         assert set(result[0]) == {1}
         assert all(a[i].all() for i in result)
 
+        # make sure wrapping works if all values are True
+        arr = g.np.ones(10, dtype=bool)
+        result = blocks(arr, min_len=1, wrap=True, only_nonzero=True)
+        assert len(result) == 1
+        assert set(result[0]) == set(range(10))
+
+        # and all false
+        arr = g.np.zeros(10, dtype=bool)
+        result = blocks(arr, min_len=1, wrap=True, only_nonzero=True)
+        assert len(result) == 0
+
+        arr = g.np.zeros(10, dtype=bool)
+        result = blocks(arr, min_len=1, wrap=True, only_nonzero=False)
+        assert len(result) == 1
+        assert set(result[0]) == set(range(10))
+
     def test_block_wrap(self):
         """
         Test blocks with wrapping
diff --git a/tests/test_interval.py b/tests/test_interval.py
index 245c1391..8e942d11 100644
--- a/tests/test_interval.py
+++ b/tests/test_interval.py
@@ -15,25 +15,42 @@ class IntervalTest(g.unittest.TestCase):
                 [[5, 15], [7, 10]],
                 [[5, 10], [10, 9]],
                 [[0, 1], [0.9, 10]],
-            ]
+                [[1000, 1001], [2000, 2001]],
+            ],
+            dtype=g.np.float64,
         )
-        tru_hit = [False, False, False, True, True, True, True]
-        tru_int = g.np.array(
-            [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [10, 20], [7, 10], [9, 10], [0.9, 1.0]]
+
+        # true intersection ranges
+        truth = g.np.array(
+            [
+                [0.0, 0.0],
+                [0.0, 0.0],
+                [0.0, 0.0],
+                [10, 20],
+                [7, 10],
+                [9, 10],
+                [0.9, 1.0],
+                [0, 0],
+            ],
+            dtype=g.np.float64,
         )
 
-        func = g.trimesh.interval.intersection
+        intersection = g.trimesh.interval.intersection
+        union = g.trimesh.interval.union
 
         # check the single- interval results
-        for ab, h, i in zip(pairs, tru_hit, tru_int):
-            r_h, r_i = func(*ab)
-            assert g.np.allclose(r_i, i)
-            assert r_h == h
+        for ab, tru in zip(pairs, truth):
+            result = intersection(*ab)
+            assert g.np.allclose(result, tru)
 
         # check the vectorized multiple interval results
-        r_h, r_i = func(pairs[:, 0, :], pairs[:, 1, :])
-        assert g.np.allclose(r_h, tru_hit)
-        assert g.np.allclose(r_i, tru_int)
+        inter = intersection(pairs[:, 0, :], pairs[:, 1, :])
+
+        assert g.np.allclose(truth, inter)
+
+        # now just run a union on these for the fun of it
+        u = union(pairs.reshape((-1, 2)))
+        assert g.np.allclose(u, [[0.0, 21.0], [1000.0, 1001.0], [2000.0, 2001.0]])
 
 
 if __name__ == "__main__":
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 33bcf927..3aac4d58 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -78,7 +78,7 @@ class MeshTests(g.unittest.TestCase):
                 # on a Path3D object
                 test = outline.paths  # NOQA
 
-            smoothed = mesh.smoothed()  # NOQA
+            smoothed = mesh.smooth_shaded  # NOQA
 
             assert abs(mesh.volume) > 0.0
 
diff --git a/tests/test_segments.py b/tests/test_segments.py
index 45daa87b..7a023155 100644
--- a/tests/test_segments.py
+++ b/tests/test_segments.py
@@ -113,6 +113,31 @@ class SegmentsTest(g.unittest.TestCase):
         # make sure overall length hasn't changed
         assert g.np.isclose(length(res), length(seg))
 
+    def test_clean(self):
+        from trimesh.path.segments import clean, resample
+
+        seg = g.np.array(
+            [[[0, 0], [1, 0]], [[1, 0], [1, 1]], [[1, 1], [0, 1]], [[0, 1], [0, 0]]],
+            dtype=g.np.float64,
+        )
+
+        c = clean(seg)
+        assert len(seg) == len(c)
+        # bounding box should be the same
+        assert g.np.allclose(
+            c.reshape((-1, 2)).min(axis=0), seg.reshape((-1, 2)).min(axis=0)
+        )
+        assert g.np.allclose(
+            c.reshape((-1, 2)).max(axis=0), seg.reshape((-1, 2)).max(axis=0)
+        )
+
+        # resample to shorten
+        r = resample(seg, maxlen=0.3)
+        assert r.shape == (16, 2, 2)
+        # after cleaning should be back to 4 segments
+        rc = clean(r)
+        assert rc.shape == (4, 2, 2)
+
     def test_svg(self):
         from trimesh.path.segments import to_svg
 
diff --git a/tests/test_smooth.py b/tests/test_smooth.py
index 6b0144e5..2f81e138 100644
--- a/tests/test_smooth.py
+++ b/tests/test_smooth.py
@@ -7,7 +7,7 @@ except BaseException:
 class SmoothTest(g.unittest.TestCase):
     def test_smooth(self):
         m = g.get_mesh("chair_model.obj", force="mesh")
-        s = m.smoothed()
+        s = m.smooth_shaded
 
         ori = g.np.hstack((m.visual.uv, m.vertices))
         check = g.np.hstack((s.visual.uv, s.vertices))