Repository: mikedh/trimesh
Pull request #2006 (2023-08-18 14:53:54)

Preserve texture on slice_plane
For the `slice_plane` method, I interpolate the UVs for the new vertices by retrieving the barycentric coordinates of the point wrt. the original face.

I'm not sure what to do when `cap=True`. For now, I'm just ignoring it and not returning any UVs when this happens.

An example:

<img width="500" alt="Screenshot 2023-08-18 at 16 34 53" src="https://github.com/mikedh/trimesh/assets/9162284/b1a16488-bdcb-446c-9d87-e3bc9ce81b54">
<img width="500" alt="Screenshot 2023-08-18 at 16 35 43" src="https://github.com/mikedh/trimesh/assets/9162284/48c07008-2d3f-45ce-bc92-d3b8d159c412">



Patch:
diff --git a/Dockerfile b/Dockerfile
index 3f1e259b..4a47ab7a 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -5,10 +5,7 @@ LABEL maintainer="mikedh@kerfed.com"
 COPY --chmod=755 docker/trimesh-setup /usr/local/bin/
 
 # Create a local non-root user.
-RUN useradd -m \ # -m creates a home directory
-    -u 499 \ # set UID to that of a system user 100-499
-    -s /bin/bash \ # set default shell to bash
-    user # make the username `user`
+RUN useradd -m -s /bin/bash user
 
 # Required for Python to be able to find libembree.
 ENV LD_LIBRARY_PATH="/usr/local/lib:$LD_LIBRARY_PATH"
@@ -30,11 +27,7 @@ COPY --chown=user:user setup.py /home/user/
 USER user
 
 # install trimesh into .local
-# then delete any included tests
-# also remove Cython after all the building is complete
-RUN pip install /home/user[all] && \
-    find /home/user/.local -type d -name tests -prune -exec rm -rf {} \; && \
-    pip uninstall -y cython
+RUN pip install /home/user[all]
 
 ####################################
 ### Build output image most things should run on
diff --git a/setup.py b/setup.py
index b8ad1795..68f72d90 100644
--- a/setup.py
+++ b/setup.py
@@ -38,6 +38,7 @@ requirements_easy = set([
     'shapely',   # handle 2D polygons robustly
     'rtree',     # create N-dimension trees for broad-phase queries
     'svg.path',  # handle SVG format path strings
+    'sympy',     # do analytical math
     'pillow',    # load images
     'embreex',   # Intel's Embree ray check engine with wheels
     'requests',  # do network requests
@@ -55,6 +56,8 @@ requirements_easy = set([
 requirements_all = requirements_easy.union([
     'python-fcl',    # do fast 3D collision queries
     'psutil',        # figure out how much memory we have
+    'glooey',        # make GUI applications with 3D stuff
+    'meshio',        # load a number of additional mesh formats; Python 3.5+
     'scikit-image',  # marching cubes and other nice stuff
     'xatlas',        # texture unwrapping
     'pyglet<2',  # render preview windows nicely : note pyglet 2.0 is basically a re-write
diff --git a/trimesh/integrate.py b/trimesh/integrate.py
new file mode 100644
index 00000000..65cd7d11
--- /dev/null
+++ b/trimesh/integrate.py
@@ -0,0 +1,129 @@
+"""
+integrate.py
+--------------
+
+Utilities for integrating functions over meshes surfaces.
+"""
+
+import sympy as sp
+from sympy.parsing.sympy_parser import parse_expr as sympy_parse
+
+from .constants import log
+from . import util
+
+import warnings
+warnings.warn(
+    '`trimesh.integrate` is deprecated and will ' +
+    'be removed in March 2024, you might take a ' +
+    'look at `quadpy` which is far more complete',
+    category=DeprecationWarning,
+    stacklevel=2)
+
+
+def symbolic_barycentric(function):
+    """
+    Symbolically integrate a function(x,y,z) across a triangle or mesh.
+
+    Parameters
+    ----------
+    function: string or sympy expression
+              x, y, z will be replaced with a barycentric representation
+              and the the function is integrated across the triangle.
+
+    Returns
+    ----------
+    evaluator: numpy lambda function of result which takes a mesh
+    expr:      sympy expression of result
+
+    Examples
+    -----------
+
+    In [1]: function = '1'
+
+    In [2]: integrator, expr = integrate_barycentric(function)
+
+    In [3]: integrator
+    Out[3]: <__main__.evaluator instance at 0x7f66cd2a6200>
+
+    In [4]: expr
+    Out[4]: 1/2
+
+    In [5]: result  = integrator(mesh)
+
+    In [6]: mesh.area
+    Out[6]: 34.641016151377542
+
+    In [7]: result.sum()
+    Out[7]: 34.641016151377542
+    """
+
+    class evaluator:
+
+        def __init__(self, expr, expr_args):
+            self.lambdified = sp.lambdify(args=expr_args,
+                                          expr=expr,
+                                          modules=['numpy', 'sympy'])
+
+        def __call__(self, mesh, *args):
+            """
+            Quickly evaluate the surface integral across a mesh
+
+            Parameters
+            ----------
+            mesh: Trimesh object
+
+            Returns
+            ----------
+            integrated: (len(faces),) float, integral evaluated for each face
+            """
+            integrated = self.lambdified(*mesh.triangles.reshape((-1, 9)).T)
+            integrated *= 2 * mesh.area_faces
+            return integrated
+
+    function, symbols = substitute_barycentric(function)
+
+    b1, b2, x1, x2, x3, y1, y2, y3, z1, z2, z3 = symbols
+    # do the first integral for b1
+    integrated_1 = sp.integrate(function, b1)
+    integrated_1 = (integrated_1.subs({b1: 1 - b2}) -
+                    integrated_1.subs({b1: 0}))
+
+    integrated_2 = sp.integrate(integrated_1, b2)
+    integrated_2 = (integrated_2.subs({b2: 1}) -
+                    integrated_2.subs({b2: 0}))
+
+    lambdified = evaluator(expr=integrated_2,
+                           expr_args=[x1, y1, z1, x2, y2, z2, x3, y3, z3])
+
+    return lambdified, integrated_2
+
+
+def substitute_barycentric(function):
+    if util.is_string(function):
+        function = sympy_parse(function)
+    # barycentric coordinates
+    b1, b2 = sp.symbols('b1 b2',
+                        real=True,
+                        positive=True)
+    # vertices of the triangles
+    x1, x2, x3, y1, y2, y3, z1, z2, z3 = sp.symbols(
+        'x1,x2,x3,y1,y2,y3,z1,z2,z3', real=True)
+
+    # generate the substitution dictionary to convert from cartesian to barycentric
+    # since the input could have been a sympy expression or a string
+    # that we parsed substitute based on name to avoid id(x) issues
+    substitutions = {}
+    for symbol in function.free_symbols:
+        if symbol.name == 'x':
+            substitutions[symbol] = b1 * x1 + b2 * x2 + (1 - b1 - b2) * x3
+        elif symbol.name == 'y':
+            substitutions[symbol] = b1 * y1 + b2 * y2 + (1 - b1 - b2) * y3
+        elif symbol.name == 'z':
+            substitutions[symbol] = b1 * z1 + b2 * z2 + (1 - b1 - b2) * z3
+    # apply the conversion to barycentric
+    function = function.subs(substitutions)
+    log.debug('converted function to barycentric: %s', str(function))
+
+    symbols = (b1, b2, x1, x2, x3, y1, y2, y3, z1, z2, z3)
+
+    return function, symbols
diff --git a/trimesh/intersections.py b/trimesh/intersections.py
index 8dd2f771..2784e50c 100644
--- a/trimesh/intersections.py
+++ b/trimesh/intersections.py
@@ -13,6 +13,7 @@ from . import triangles as tm
 from . import transformations as tf
 
 from .constants import tol
+from .triangles import points_to_barycentric
 
 
 def mesh_plane(mesh,
@@ -430,6 +431,7 @@ def slice_faces_plane(vertices,
                       faces,
                       plane_normal,
                       plane_origin,
+                      uv=None,
                       face_index=None,
                       cached_dots=None):
     """
@@ -447,6 +449,8 @@ def slice_faces_plane(vertices,
         Normal vector of plane to intersect with mesh
     plane_origin :  (3,) float
         Point on plane to intersect with mesh
+    uv : (n, 2) float, optional
+        UV coordinates of source mesh to slice
     face_index : ((m,) int)
         Indexes of faces to slice. When no mask is provided, the
         default is to slice all faces.
@@ -460,10 +464,14 @@ def slice_faces_plane(vertices,
         Vertices of sliced mesh
     new_faces : (n, 3) int
         Faces of sliced mesh
+    new_uv : (n, 2) int or None
+        UV coordinates of sliced mesh
     """
 
     if len(vertices) == 0:
-        return vertices, faces
+        return vertices, faces, uv
+
+    have_uv = uv is not None
 
     # Construct a mask for the faces to slice.
     if face_index is not None:
@@ -539,7 +547,8 @@ def slice_faces_plane(vertices,
         if len(new_faces) == 0:
             # if no new faces at all return empty arrays
             empty = (np.zeros((0, 3), dtype=np.float64),
-                     np.zeros((0, 3), dtype=np.int64))
+                     np.zeros((0, 3), dtype=np.int64),
+                     np.zeros((0, 2), dtype=np.float64) if have_uv else None)
             return empty
 
         # find the unique indices in the new faces
@@ -551,8 +560,9 @@ def slice_faces_plane(vertices,
         # use the unique indices for our final vertices and faces
         final_vert = vertices[unique]
         final_face = inverse.reshape((-1, 3))
+        final_uv = uv[unique] if have_uv else None
 
-        return final_vert, final_face
+        return final_vert, final_face, final_uv
 
     # Extract the intersections of each triangle's edges with the plane
     o = cut_triangles                               # origins
@@ -566,6 +576,8 @@ def slice_faces_plane(vertices,
 
     # Initialize the array of new vertices with the current vertices
     new_vertices = vertices
+    new_quad_vertices = np.zeros((0, 3))
+    new_tri_vertices = np.zeros((0, 3))
 
     # Handle the case where a new quad is formed by the intersection
     # First, extract the intersection points belonging to a new quad
@@ -641,7 +653,23 @@ def slice_faces_plane(vertices,
     final_vert = new_vertices[unique]
     final_face = inverse.reshape((-1, 3))
 
-    return final_vert, final_face
+    final_uv = None
+    if have_uv:
+        # Generate barycentric coordinates for intersection vertices
+        quad_barycentrics = points_to_barycentric(
+            np.repeat(vertices[cut_faces_quad], 2, axis=0),
+            new_quad_vertices)
+        tri_barycentrics = points_to_barycentric(
+            np.repeat(vertices[cut_faces_tri], 2, axis=0),
+            new_tri_vertices)
+        all_barycentrics = np.concatenate([quad_barycentrics, tri_barycentrics])
+
+        # Interpolate UVs
+        cut_uv = np.concatenate([uv[cut_faces_quad], uv[cut_faces_tri]])
+        new_uv = np.einsum('ijk,ij->ik', np.repeat(cut_uv, 2, axis=0), all_barycentrics)
+        final_uv = np.concatenate([uv, new_uv])[unique]
+
+    return final_vert, final_face, final_uv
 
 
 def slice_mesh_plane(mesh,
@@ -689,6 +717,7 @@ def slice_mesh_plane(mesh,
 
     # avoid circular import
     from .base import Trimesh
+    from .visual import TextureVisuals
     from .path import polygons
     from scipy.spatial import cKDTree
     from .creation import triangulate_polygon
@@ -712,6 +741,11 @@ def slice_mesh_plane(mesh,
     vertices = mesh.vertices.copy()
     faces = mesh.faces.copy()
 
+    # We copy the UV coordinates if available
+    has_uv = (hasattr(mesh.visual, 'uv') and np.shape(
+        mesh.visual.uv) == (len(mesh.vertices), 2)) and not cap
+    uv = mesh.visual.uv.copy() if has_uv else None
+
     if 'process' not in kwargs:
         kwargs['process'] = False
 
@@ -719,9 +753,10 @@ def slice_mesh_plane(mesh,
     for origin, normal in zip(plane_origin.reshape((-1, 3)),
                               plane_normal.reshape((-1, 3))):
         # save the new vertices and faces
-        vertices, faces = slice_faces_plane(
+        vertices, faces, uv = slice_faces_plane(
             vertices=vertices,
             faces=faces,
+            uv=uv,
             plane_normal=normal,
             plane_origin=origin,
             face_index=face_index)
@@ -774,5 +809,8 @@ def slice_mesh_plane(mesh,
                 faces.append(vid[fn])
             faces = np.vstack(faces)
 
+    visual = TextureVisuals(
+        uv=uv, material=mesh.visual.material.copy()) if has_uv else None
+
     # return the sliced mesh
-    return Trimesh(vertices=vertices, faces=faces, **kwargs)
+    return Trimesh(vertices=vertices, faces=faces, visual=visual, **kwargs)


Test patch:
diff --git a/tests/test_section.py b/tests/test_section.py
index 2f804497..a75ea950 100644
--- a/tests/test_section.py
+++ b/tests/test_section.py
@@ -291,6 +291,32 @@ class SliceTest(g.unittest.TestCase):
 
         assert len(sliced.faces) > 0
 
+    def test_textured_mesh(self):
+        # Create a plane mesh with UV == xy
+        plane = g.trimesh.Trimesh(
+            vertices=g.np.array([[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.]]),
+            faces=g.np.array([[0, 1, 2], [2, 3, 0]]))
+        plane.visual = g.trimesh.visual.TextureVisuals(
+            uv=plane.vertices[:, :2], material=g.trimesh.visual.material.empty_material())
+
+        # We cut the plane and check that the new UV match the new vertices
+        origin = g.np.array([0.5, 0.5, 0.])
+        normal = g.trimesh.unitize([2, 1, 2])
+        sliced = plane.slice_plane(plane_origin=origin, plane_normal=normal)
+        assert g.np.isclose(sliced.vertices[:, :2], sliced.visual.uv).all()
+
+        # Test scenario when plane does not cut
+        origin = g.np.array([-1., -1., 0.])
+        normal = g.trimesh.unitize([1, 1, 2])
+        sliced = plane.slice_plane(plane_origin=origin, plane_normal=normal)
+        assert g.np.isclose(sliced.vertices[:, :2], sliced.visual.uv).all()
+
+        # Test cut with no new vertices
+        origin = g.np.array([0.5, 0.5, 0.])
+        normal = g.trimesh.unitize([2, -2, 1])
+        sliced = plane.slice_plane(plane_origin=origin, plane_normal=normal)
+        assert g.np.isclose(sliced.vertices[:, :2], sliced.visual.uv).all()
+
     def test_cap_coplanar(self):
         # check to see if we handle capping with
         # existing coplanar faces correctly