Performance-focused rewrite of the ConvexPolyhedron class

[janbridley]

Prior to this PR, the ConvexPolyhedron class used Scipy's (Qhull) convex hull to compute faces of the solid. The shape was then initialized in the same manner as a standard Polyhedron. This rewrite stores simplices and other data from the convex hull, which can be used to accelerate many of the classes methods. Slow and commonly run functions have been rewritten for performance, and new and/or improved pytests have been written.

Description

The new class builds a convex hull, then stores most of the data from that object. This data, especially the list of simplices, allow for numpy vectorization of many functions that required a for loop in the previous code. These simplices can be sorted to ensure math involving the curl/divergence theorems works properly, and can be merged based on common normals to build faces similar to the current merge_faces function.

The following methods have been reimplemented, and a brief description is given:

• find_equations: vectorized using numpy
• _surface_triangulation: updated to use precalculated simplices, so polytri is not required
• volume, surface_area: reimplemented to use stored volume data for the calculations.
• centroid: rewritten to calculate from a list of simplices. Vectorized using numpy, and computes solutions using (_centroid_from_triangulated_surface)
• _compute_inertia_tensor: fully rewritten for performance, using an algorithm that works with tetrahedral decompositions that include negative-volume tetrahedra. It is also significantly faster.
• _rescale: updated to properly scale all stored properties. Also recalculates the centroid to ensure proper behavior.
• get_face_area: vectorized with numpy

The following new methods have been added:

• _combine simplices: faster alternative to merge faces, but only works on simplices
• _find_simplex_equations: similar to _find_equations, but works on every simplex rather than every face
• _find_coplanar_simplices: helper function that finds simplices with the same plane equation (within a tolerance).
• _calculate_signed_volume: faster method of computing volume via a tetrahedral decomposition. Used as an internal check that volume is stable when modifying the shape
• _find_face_centroids: used with the new face_centroids property. Quickly finds the center of mass for each of the polygonal faces

The following properties have been added:

• simplices
• equations
• face_centroids

Motivation and Context

Creation of a ConvexPolyhedron object takes on the order of 1e-2s - this is fine when generating small numbers of shapes, but can become a significant hindrance when generating tens or hundreds of thousands of shapes. In addition, many commonly used methods (especially volume, centroid, and inertia tensor) took on the order of 1e-3s and did not cache their data. These changes address these performance issues, and add a wide range of private variables which can be used to more quickly calculate geometric properties that have not yet been implemented. Benchmarks are included below, but in general this rewrite provides a 10x-100x performance improvement on improved methods. For properties with caching, values are only recalculated when the shape is modified and access the data takes on the order of 1e-8s

Benchmarks

The chart above shows the runtimes of common actions averages over 10 thousand samples of each Johnson solid. Errorbars show the standard deviation between shapes.

Note that tests were performed in much more depth than shown below - these are simply a subset of the total

With optimization

Cube (V=8, F=6, E=12)

cube = PlatonicFamily.get_shape("Cube")

%timeit -r 10 -n 10_000 ConvexPolyhedron(cube.vertices)
> 301 us ± 34.2 us per shape

%timeit -r 10 -n 100_000 cube._calculate_signed_volume()
> 7.69 us ± 239 ns per shape

%timeit -r 10 -n 100_000 cube._centroid_from_triangulated_surface()
> 25.7 us ± 676 ns per shape

%timeit -r 10 -n 100_000 cube.inertia_tensor
> 111 us ± 2.74 us per shape

Metabidiminished rhombicosidodecahedron (V=50, F=42, E=90)

mbri = JohnsonFamily.get_shape("Metabidiminished Rhombicosidodecahedron")

%timeit -r 10 -n 10_000 ConvexPolyhedron(mbri.vertices)
> 1.81 ms ± 315 us per shape

%timeit -r 10 -n 100_000 mbri._calculate_signed_volume()
> 20.6 us ± 2.38 us per shape

%timeit -r 10 -n 100_000 mbri._centroid_from_triangulated_surface()
> 32.2 us ± 1.96 us per shape

%timeit -r 10 -n 100_000 mbri.inertia_tensor
> 194 us ± 9.65 us per shape

Without optimization (previous versions of coxeter)

Cube (V=8, F=6, E=12)

cube = PlatonicFamily.get_shape("Cube")

%timeit -r 10 -n 10_000 ConvexPolyhedron(cube.vertices)
> 5.22 ms ± 169 us per shape

%timeit -r 10 -n 100_000 cube.volume
> 1.38 ms ± 22.8 us per shape

%timeit -r 10 -n 100_000 cube.centroid
> 946 us ± 37.4 us per shape

%timeit -r 10 -n 100_000 cube.inertia_tensor
> 5.32 ms ± 478 us per shape

Metabidiminished rhombicosidodecahedron (V=50, F=42, E=90)

mbri = JohnsonFamily.get_shape("Metabidiminished Rhombicosidodecahedron")

%timeit -r 10 -n 1_0000 ConvexPolyhedron(mbri.vertices)
> 38.1 ms ± 1.05 ms per shape

%timeit -r 10 -n 10_000 mbri.volume
> 10.2 ms ± 328 us per shape

%timeit -r 10 -n 10_000 mbri.centroid
> 8.34 ms ± 829 us per shape

%timeit -r 10 -n 10_000 mbri.inertia_tensor
> 34.5 ms ± 1 ms per shape

Types of Changes

• Documentation update
• Bug fix
• New feature
• Breaking change¹

¹The change breaks (or has the potential to break) existing functionality.

Checklist:

• I am familiar with the Contributing Guidelines.
• I agree with the terms of the Contributor Agreement.
• The changes introduced by this pull request are covered by existing or newly introduced tests.
• I have updated the changelog and the credits.

[DomFijan]

Is this necessary? This just changes which vertex you start at basically
right? I might be missing something. arctan2 codomain is -pi to pi.

[janbridley]

Pytests for ConvexPolyhedron pass with or without the shift, but the current code mirrors the behavior of shapes/_reorder_verts. Removing the rescaling from ConvexPolygon causes the test test_polygon.p::test_nonregular_convex_polygon_distance_to_surface_unit_area_ngon to fail. It seems that reordering the vertices is somehow required to ensure the resulting polygon is simple. On a first glance I can't figure out why this would be the case, but I think we should explore further because the behavior is quite odd. This also highlights that we should add some 3D test cases with non regular faces - the current named_solids_mark does not include anything more irregular than Johnson solids

While the above comment is mostly correct, the damasceno_shapes_mark does test seven irregular polyhedra (O15-O22). The EllipsoidSurfaceStrategy also tests irregular shapes, although almost all of them will be deltahedra

[DomFijan]

I've left some comments that we should discuss. Some methods are lacking tests that should be added. And some minor tweaks are suggested.

[DomFijan]

I don't have any more suggestions for improvement.This looks fine and since tests pass I think this might be good to go.

[janbridley]

@tcmoore3 Any other feedback? If not I will merge in the morning and start working on a release