-
Notifications
You must be signed in to change notification settings - Fork 633
Description
Hello,
We ran into an edge case in our software and tracked it down to the mesh.projected(precise=True) call. Now this is a long story, so I'm going to cut to the chase. We save our meshes as VTP files using PyVista, then read them later using PyVista and convert to a Trimesh object using trimesh.Trimesh(pv_mesh.points, pv_mesh.regular_faces). After this, we call mesh.projected with a set of very specific normals. With a very specific single mesh at a specific single normal, this fails. We have tried many thousands of other meshes with these normals before without any problems, so this is extremely niche. The mesh that fails is uploaded in the example below (54 vertices and 104 faces). The normal that fails with this mesh is also in the example below.
I am not sure if this a PyVista problem with the vertex positions array (see workarounds), if this is a problem with Trimesh, or a problem with shapely. I decided to start here and we can move the issue if desired.
Minimal Example
feature.zip (download and extract the single file it contains, GitHub won't let me upload a VTP file)
import trimesh
import pyvista as pv
pv_mesh: pv.PolyData = pv.read("feature.vtp")
mesh = trimesh.Trimesh(pv_mesh.points, pv_mesh.regular_faces)
mesh.projected([-0.7941530735495742, 0.5769859816799351, 0.19080899537654492], precise=True) # CrashesGEOSException: TopologyException: Ring edge missing at -1.5732562941774 7.904434391775184
Stack trace
---------------------------------------------------------------------------
GEOSException Traceback (most recent call last)
...
File .venv/lib/python3.13/site-packages/trimesh/base.py:2726, in Trimesh.projected(self, normal, **kwargs)
2723 from .path import Path2D
2724 from .path.polygons import projected
-> 2726 projection = projected(mesh=self, normal=normal, **kwargs)
2727 if projection is None:
2728 return Path2D()
File .venv/lib/python3.13/site-packages/trimesh/path/polygons.py:836, in projected(mesh, normal, origin, ignore_sign, rpad, apad, tol_dot, precise)
833 faces = mesh.faces[side]
834 # just union all the polygons
835 return (
--> 836 ops.unary_union(
837 [Polygon(f) for f in vertices_2D[np.column_stack((faces, faces[:, :1]))]]
838 )
839 .buffer(eps)
840 .buffer(-eps)
841 )
843 # a sequence of face indexes that are connected
844 face_groups = graph.connected_components(adjacency, nodes=np.nonzero(side)[0])
File .venv/lib/python3.13/site-packages/shapely/ops.py:117, in CollectionOperator.unary_union(self, geoms)
111 def unary_union(self, geoms):
112 """Return the union of a sequence of geometries.
113
114 Usually used to convert a collection into the smallest set of polygons
115 that cover the same area.
116 """
--> 117 return shapely.union_all(geoms, axis=None)
File .venv/lib/python3.13/site-packages/shapely/decorators.py:173, in deprecate_positional.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
171 @wraps(func)
172 def wrapper(*args, **kwargs):
--> 173 result = func(*args, **kwargs)
175 n = len(args)
176 if n > warn_from:
File .venv/lib/python3.13/site-packages/shapely/decorators.py:88, in multithreading_enabled.<locals>.wrapped(*args, **kwargs)
86 for arr in array_args:
87 arr.flags.writeable = False
---> 88 return func(*args, **kwargs)
89 finally:
90 for arr, old_flag in zip(array_args, old_flags):
File .venv/lib/python3.13/site-packages/shapely/set_operations.py:553, in union_all(geometries, grid_size, axis, **kwargs)
549 raise ValueError("grid_size parameter only accepts scalar values")
551 return lib.unary_union_prec(collections, grid_size, **kwargs)
--> 553 return lib.unary_union(collections, **kwargs)
GEOSException: TopologyException: Ring edge missing at -1.5732562941774 7.904434391775184
Workarounds
- Change
pv_mesh.pointstopv_mesh.points.astype(np.float32) - Change
pv_mesh.pointstonp.round(pv_mesh.points, 64) - Saving as STL and load that instead