From f0b8598e445f50fdbf9bbdfc887b9b7c6cb4a08d Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Wed, 25 Jan 2023 18:30:33 -0500 Subject: [PATCH 1/4] test: Fix assumption on quaternions. Only test quaternions that are not effectively zero. --- tests/test_polygon.py | 4 ++-- tests/test_spheropolygon.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_polygon.py b/tests/test_polygon.py index 4c6d78be..9973f785 100644 --- a/tests/test_polygon.py +++ b/tests/test_polygon.py @@ -189,7 +189,7 @@ def test_convex_area(points): @given(random_quat=arrays(np.float64, (4,), elements=floats(-1, 1, width=64))) def test_rotation_signed_area(random_quat): """Ensure that rotating does not change the signed area.""" - assume(not np.all(random_quat == 0)) + assume(not np.allclose(random_quat, 0)) random_quat = rowan.normalize(random_quat) rotated_points = rowan.rotate(random_quat, get_square_points()) poly = Polygon(rotated_points) @@ -253,7 +253,7 @@ def test_bounding_circle_radius_random_hull(points): ) def test_bounding_circle_radius_random_hull_rotation(points, rotation): """Test that rotating vertices does not change the bounding radius.""" - assume(not np.all(rotation == 0)) + assume(not np.allclose(rotation, 0)) hull = ConvexHull(points) poly = Polygon(points[hull.vertices]) diff --git a/tests/test_spheropolygon.py b/tests/test_spheropolygon.py index a9df4a0e..a1d34b6e 100644 --- a/tests/test_spheropolygon.py +++ b/tests/test_spheropolygon.py @@ -152,7 +152,7 @@ def test_convex_signed_area(square_points): ) ) def testfun(random_quat): - assume(not np.all(random_quat == 0)) + assume(not np.allclose(random_quat, 0)) random_quat = rowan.normalize(random_quat) rotated_points = rowan.rotate(random_quat, square_points) r = 1 From 9e6c0885284e37796e8bd42b68b12a79dae1a451 Mon Sep 17 00:00:00 2001 From: janbridley Date: Mon, 20 Feb 2023 21:03:53 -0500 Subject: [PATCH 2/4] Changed miniball solver to a faster and more accurate variant Replacing the ``miniball`` package with the improved ``cyminiball`` c++ binding allows for accurate calculation of a minimal_bounding_sphere for complex polyhedra. An internal test to ensure vertices are within (or extremely close to) the minimal bounding sphere confirms expected behavior and catches cases where an incorrect miniball could be calculated. --- coxeter/shapes/polygon.py | 10 +++++----- coxeter/shapes/polyhedron.py | 31 ++++++++++++++++++++++--------- setup.py | 2 +- 3 files changed, 28 insertions(+), 15 deletions(-) diff --git a/coxeter/shapes/polygon.py b/coxeter/shapes/polygon.py index 6edbdda3..c5be5c7b 100644 --- a/coxeter/shapes/polygon.py +++ b/coxeter/shapes/polygon.py @@ -16,7 +16,7 @@ from .utils import _generate_ax, rotate_order2_tensor, translate_inertia_tensor try: - import miniball + import cyminiball MINIBALL = True except ImportError: @@ -475,13 +475,13 @@ def minimal_bounding_circle(self): """:class:`~.Circle`: Get the minimal bounding circle.""" if not MINIBALL: raise ImportError( - "The miniball module must be installed. It can " + "The cyminiball module must be installed. It can " "be installed as an extra with coxeter (e.g. " "with pip install coxeter[bounding_sphere], or " - "directly from PyPI using pip install miniball." + "directly from PyPI using pip install cyminiball." ) - # The algorithm in miniball involves solving a linear system and + # The algorithm in cyminiball involves solving a linear system and # can therefore occasionally be somewhat unstable. Applying a # random rotation will usually fix the issue. max_attempts = 10 @@ -491,7 +491,7 @@ def minimal_bounding_circle(self): while attempt < max_attempts: attempt += 1 try: - center, r2 = miniball.get_bounding_ball(vertices) + center, r2 = cyminiball.compute(vertices) break except np.linalg.LinAlgError: current_rotation = rowan.random.rand(1) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 074db4b6..da4d02c8 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -18,7 +18,7 @@ from .utils import _generate_ax, _set_3d_axes_equal, translate_inertia_tensor try: - import miniball + import cyminiball MINIBALL = True except ImportError: @@ -550,27 +550,40 @@ def minimal_bounding_sphere(self): """:class:`~.Sphere`: Get the polyhedron's bounding sphere.""" if not MINIBALL: raise ImportError( - "The miniball module must be installed. It can " + "The cyminiball module must be installed. It can " "be installed as an extra with coxeter (e.g. " 'with "pip install coxeter[bounding_sphere]") or ' - 'directly from PyPI using "pip install miniball".' + 'directly from PyPI using "pip install cyminiball".' ) - - # The algorithm in miniball involves solving a linear system and - # can therefore occasionally be somewhat unstable. Applying a - # random rotation will usually fix the issue. - max_attempts = 10 + # The miniball algorithm (in this case, Gärtner's miniball) can + # experience numerical instability of ~10x machine epsilon. If an incorrect + # miniball is found, applying random rotations will usually fix the issue. + + # The following subfunction checks if all vertices lie within the + # minimal bounding sphere, to a tolerance of 1e-15. + def verify_inside(points, center, radius, eps=1e-15): + return np.all(np.linalg.norm(points - center, axis=1) <= radius + eps) + + # Simple polyhedra (e.g. dodecahedra) will require far less than 200 attempts + # to compute a correct miniball. Polyhedra with augmentations and large numbers + # of vertices will take more, but none are likely to exceed 200 attempts. Worst- + # case runtime is approximately 0.25s, and occurrs in ~5/1e5 samples. + max_attempts = 200 attempt = 0 current_rotation = [1, 0, 0, 0] vertices = self.vertices while attempt < max_attempts: attempt += 1 try: - center, r2 = miniball.get_bounding_ball(vertices) + center, r2 = cyminiball.compute(vertices) + assert verify_inside(vertices, center, np.sqrt(r2)) break except np.linalg.LinAlgError: current_rotation = rowan.random.rand(1) vertices = rowan.rotate(current_rotation, vertices) + except AssertionError: + current_rotation = rowan.random.rand(1) + vertices = rowan.rotate(current_rotation, vertices) else: raise RuntimeError("Unable to solve for a bounding sphere.") diff --git a/setup.py b/setup.py index 22a45faa..891bc8bd 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ ] bounding_deps = [ - "miniball", + "cyminiball", ] extras = { From df17105311e9d0355270433b2ee27b9a23fcb466 Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Tue, 21 Feb 2023 08:39:44 -0500 Subject: [PATCH 3/4] Fixed typo. Co-authored-by: Bradley Dice --- coxeter/shapes/polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index da4d02c8..e9955022 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -567,7 +567,7 @@ def verify_inside(points, center, radius, eps=1e-15): # Simple polyhedra (e.g. dodecahedra) will require far less than 200 attempts # to compute a correct miniball. Polyhedra with augmentations and large numbers # of vertices will take more, but none are likely to exceed 200 attempts. Worst- - # case runtime is approximately 0.25s, and occurrs in ~5/1e5 samples. + # case runtime is approximately 0.25s, and occurs in ~5/1e5 samples. max_attempts = 200 attempt = 0 current_rotation = [1, 0, 0, 0] From a69457c17dcc4fb7e525752218e7f83f6f4dcaa5 Mon Sep 17 00:00:00 2001 From: Jen Bradley <55467578+janbridley@users.noreply.github.com> Date: Tue, 21 Feb 2023 08:41:07 -0500 Subject: [PATCH 4/4] Changed miniball random rotations to apply to original vertices Co-authored-by: Bradley Dice --- coxeter/shapes/polyhedron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index e9955022..5f79cfe6 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -583,7 +583,7 @@ def verify_inside(points, center, radius, eps=1e-15): vertices = rowan.rotate(current_rotation, vertices) except AssertionError: current_rotation = rowan.random.rand(1) - vertices = rowan.rotate(current_rotation, vertices) + vertices = rowan.rotate(current_rotation, self.vertices) else: raise RuntimeError("Unable to solve for a bounding sphere.")