nsphere speedup

tests/test_bounds.py

@@ -120,8 +120,13 @@ def test_2D(self):
     def test_cylinder(self):
         """
         """
+        # not rotationally symmetric
+        mesh = g.get_mesh('featuretype.STL')
+
         height = 10.0
         radius = 1.0
+
+        # spherical coordinates to loop through
         sphere = g.trimesh.util.grid_linspace([[0, 0],
                                                [g.np.pi * 2, g.np.pi * 2]],
                                               5)
@@ -137,6 +142,14 @@ def test_cylinder(self):
                                 p.bounding_cylinder.primitive.height,
                                 rtol=.01)
 
+            # regular mesh should have the same bounding cylinder
+            # regardless of transform
+            copied = mesh.copy()
+            copied.apply_transform(T)
+            assert g.np.isclose(mesh.bounding_cylinder.volume,
+                                copied.bounding_cylinder.volume,
+                                rtol=.05)
+
     def test_obb_order(self):
         # make sure our sorting and transform flipping of
         # OBB extents are working by checking against a box

trimesh/bounds.py

@@ -219,7 +219,7 @@ def oriented_bounds(obj, angle_digits=1, ordered=True):
     return to_origin, min_extents
 
 
-def minimum_cylinder(obj, sample_count=10, angle_tol=.001):
+def minimum_cylinder(obj, sample_count=6, angle_tol=.001):
     """
     Find the approximate minimum volume cylinder which contains
     a mesh or a a list of points.
@@ -277,7 +277,7 @@ def volume_from_angles(spherical, return_data=False):
         height = projected[:, 2].ptp()
 
         try:
-            center_2D, radius = nsphere.minimum_nsphere(projected[:, 0:2])
+            center_2D, radius = nsphere.minimum_nsphere(projected[:, :2])
         except BaseException:
             # in degenerate cases return as infinite volume
             return np.inf
@@ -297,11 +297,17 @@ def volume_from_angles(spherical, return_data=False):
 
     # sample a hemisphere so local hill climbing can do its thing
     samples = util.grid_linspace([[0, 0], [np.pi, np.pi]], sample_count)
-    # add the principal inertia vectors if we have a mesh
+
+    # if it's rotationally symmetric the bounding cylinder
+    # is almost certainly along one of the PCI vectors
+    # if hasattr(obj, 'symmetry_axis') and obj.symmetry_axis is not None:
+    #    samples = util.vector_to_spherical(obj.principal_inertia_vectors)
     if hasattr(obj, 'principal_inertia_vectors'):
+        # add the principal inertia vectors if we have a mesh
         samples = np.vstack(
-            (samples, util.vector_to_spherical(
-                obj.principal_inertia_vectors)))
+            (samples,
+             util.vector_to_spherical(obj.principal_inertia_vectors)))
+
     tic = [time.time()]
     # the projected volume at each sample
     volumes = np.array([volume_from_angles(i) for i in samples])
@@ -314,7 +320,7 @@ def volume_from_angles(spherical, return_data=False):
     step = 2 * np.pi / sample_count
     bounds = [(best[0] - step, best[0] + step),
               (best[1] - step, best[1] + step)]
-    # run the optimization
+    # run the local optimization
     r = optimize.minimize(volume_from_angles,
                           best,
                           tol=angle_tol,

trimesh/nsphere.py

@@ -8,7 +8,7 @@
 import numpy as np
 
 from . import convex
-
+from . import bounds
 from .constants import log, tol
 
 try:
@@ -70,17 +70,29 @@ def minimum_nsphere(obj):
     # find the maximum radius^2 point for each of the voronoi vertices
     # this is worst case quite expensive, but we have used quick convex
     # hull methods to reduce n for this operation
-    # we are doing comparisons on the radius^2 value so as to only do a sqrt
-    # once
-    r2 = np.array([((points - v)**2).sum(axis=1).max()
-                   for v in voronoi.vertices])
+    # we are doing comparisons on the radius squared then rooting once
+
+    if len(points) * len(voronoi.vertices) > 1e7:
+        # if we have a bajillion points loop
+        radii_2 = np.array([((points - v)**2).sum(axis=1).max()
+                            for v in voronoi.vertices])
+    else:
+        # otherwise tiling is massively faster
+        # what dimension are the points
+        dim = points.shape[1]
+        v_tile = np.tile(voronoi.vertices, len(points)).reshape((-1, dim))
+        # tile points per voronoi vertex
+        p_tile = np.tile(points, (len(voronoi.vertices), 1)).reshape((-1, dim))
+        # find the maximum radius of points for each voronoi vertex
+        radii_2 = ((p_tile - v_tile)**2).sum(axis=1).reshape(
+            (len(voronoi.vertices), -1)).max(axis=1)
+
     # we want the smallest sphere, so we take the min of the radii options
-    r2_idx = r2.argmin()
-    center_v = voronoi.vertices[r2_idx]
+    radii_idx = radii_2.argmin()
 
     # return voronoi radius and center to global scale
-    radius_v = np.sqrt(r2[r2_idx]) * points_scale
-    center_v = (center_v * points_scale) + points_origin
+    radius_v = np.sqrt(radii_2[radii_idx]) * points_scale
+    center_v = (voronoi.vertices[radii_idx] * points_scale) + points_origin
 
     if radius_v > fit_R:
         return fit_C, fit_R

trimesh/version.py

@@ -1 +1 @@
-__version__ = '2.35.26'
+__version__ = '2.35.27'