Add shapemoves subpackage based on conway

coxeter/__init__.py

@@ -18,7 +18,7 @@
 
 __version__ = "0.10.0"
 
-from . import families, io, shapes
+from . import families, io, shapemoves, shapes
 from .shape_getters import from_gsd_type_shapes
 
 __all__ = ["families", "shapes", "from_gsd_type_shapes", "io"]

coxeter/families/__init__.py

@@ -75,6 +75,7 @@
 
 from .common import (
     ArchimedeanFamily,
+    CanonicalTrapezohedronFamily,
     CatalanFamily,
     JohnsonFamily,
     PlatonicFamily,
@@ -126,6 +127,7 @@
     "UniformPrismFamily",
     "UniformPyramidFamily",
     "ArchimedeanFamily",
+    "CanonicalTrapezohedronFamily",
     "CatalanFamily",
     "JohnsonFamily",
     "PrismAntiprismFamily",

coxeter/families/common.py

@@ -17,6 +17,14 @@
 from .tabulated_shape_family import TabulatedGSDShapeFamily
 
 
+def csc(theta):
+    """Compute the cosecant of a value.
+
+    :meta private:
+    """
+    return 1 / sin(theta)
+
+
 # Allows us to monkeypatch an existing method with our choice of warning
 def _deprecated_method(func, deprecated="", replacement="", reason=""):
     @wraps(func)
@@ -32,16 +40,22 @@ def wrapper(*args, **kwargs):
 
 
 def sec(theta):
-    """Return the secant of an angle."""
+    """Return the secant of an angle.
+
+    :meta private:
+    """
     return 1 / cos(theta)
 
 
 def cot(theta):
-    """Return the cotangent of an angle."""
+    """Return the cotangent of an angle.
+
+    :meta private:
+    """
     return 1 / tan(theta)
 
 
-def _make_ngon(n, z=0, area=1, angle=0):
+def _make_ngon(n, z=0.0, area=1, angle=0):
     """Make a regular n-gon with a given area, z height, and rotation angle.
 
     The initial vertex  lies on the :math:`x` axis by default, but can be rotated by
@@ -297,6 +311,79 @@ def make_vertices(cls, n):
         return vertices
 
 
+class CanonicalTrapezohedronFamily(ShapeFamily):
+    r"""The infinite family of canonical n-trapezohedra (antiprism duals).
+
+    Formulas for vertices are derived from :cite:`Rajpoot_2015` rather than via explicit
+    canonicalization to ensure the method is deterministic and fast.
+    """
+
+    @classmethod
+    def get_shape(cls, n: int):
+        r"""Generate a canonical n-antiprism of unit volume.
+
+        Args:
+            n (int): The number of vertices of the base polygons (:math:`n \geq 3`).
+
+        Returns
+        -------
+             :class:`~.ConvexPolyhedron`: The corresponding convex polyhedron.
+        """
+        return ConvexPolyhedron(cls.make_vertices(n))
+
+    @classmethod
+    def make_vertices(cls, n):
+        r"""Generate the vertices of a uniform right n-antiprism with unit volume.
+
+        Args:
+            n (int): The number of vertices of the base polygons (:math:`n \geq 3`).
+
+        Returns
+        -------
+             :math:`(n*2 + 2, 3)` :class:`numpy.ndarray` of float:
+                 The vertices of the trapezohedron.
+        """
+        r = 0.5 * csc(np.pi / n)  # Midradius for canonical trapezohedron
+        area = r**2 * n / 2 * sin(2 * np.pi / n)  # Area of center polygons
+
+        # Height of center polygons
+        c = sqrt(4 - sec(np.pi / (2 * n)) ** 2) / (4 + 8 * cos(np.pi / n))
+
+        # Height of apexes
+        z = (
+            0.25
+            * cos(np.pi / (2 * n))
+            * cot(np.pi / (2 * n))
+            * csc(3 * np.pi / (2 * n))
+            * sqrt(4 - sec(np.pi / (2 * n)) ** 2)
+        )
+
+        top_center_polygon = _make_ngon(n, c, area, angle=np.pi / n)
+        bottom_center_polygon = _make_ngon(n, -c, area)
+
+        vertices = np.concatenate(
+            [
+                top_center_polygon,
+                bottom_center_polygon,
+                [[0, 0, z], [0, 0, -z]],
+            ]
+        )
+
+        # Compute height of (canonical) trapezohedron
+        h = 1 / 4 * csc(np.pi / (2 * n)) * sqrt(2 * cos(np.pi / n) + 1) / 2
+
+        # Compute edge lengths (required to compute the volume)
+        edge_length_ratio = 1 / (2 - 2 * cos(np.pi / n))
+        short_edge_length = np.linalg.norm(
+            top_center_polygon[0] - bottom_center_polygon[0]
+        )
+        long_edge_length = short_edge_length * edge_length_ratio
+
+        vo = 2 * n / 3 * short_edge_length * long_edge_length * h
+        vertices *= cbrt(1 / vo)  # Rescale to unit volume
+        return vertices
+
+
 PlatonicFamily = TabulatedGSDShapeFamily._from_json_file(
     os.path.join(_DATA_FOLDER, "platonic.json"),
     classname="PlatonicFamily",

coxeter/shapemoves.py

@@ -0,0 +1,110 @@
+# Copyright (c) 2015-2025 The Regents of the University of Michigan.
+# This file is from the coxeter project, released under the BSD 3-Clause License.
+
+"""Shape-move operators, as defined by geometers John Conway and George Hart."""
+
+import numpy as np
+from numpy.typing import ArrayLike
+from scipy.spatial import HalfspaceIntersection
+
+from coxeter.shapes import ConvexPolyhedron
+
+
+def _truncate(poly: ConvexPolyhedron, t: float, degrees=None, filter_unique=False):
+    # Define the distance along each edge between original vertices and new vertices
+    edge_endpoint_offsets = (poly.edge_vectors[..., None] * [t, -t]).transpose(0, 2, 1)
+    # Compute new point locations from offsets and original points
+    new_vertices = poly.vertices[poly.edges] + edge_endpoint_offsets
+
+    if degrees is not None:
+        # Compute the degree of each vertex of the polyhedron
+        vertex_degrees = np.bincount(poly.edges.ravel())
+        # Mask vertices to be truncated (vertex degree in input degree set)
+        degree_mask = np.any(vertex_degrees[:, None] == np.array(degrees), axis=1)
+        vertices_to_truncate = np.where(degree_mask)[0]
+        edge_endpoint_mask = np.isin(poly.edges, vertices_to_truncate)
+
+        # Take only vertices of correct degree and merge in original vertices
+        new_vertices = np.vstack(
+            [new_vertices[edge_endpoint_mask], poly.vertices[~degree_mask]]
+        )
+
+    # Reshape back down to an N,3 array and prune duplicates
+    return np.unique(new_vertices.reshape(-1, 3), axis=0)
+
+
+def vertex_truncate(poly: ConvexPolyhedron, t: float, degrees: ArrayLike = None):
+    """Truncate the vertices of a polyhedron.
+
+    .. important::
+
+        Conway's original definition of vertex truncation is only well-defined for
+        vertices connecting a single type of incoming edge (like the Platonic and
+        Archimedean solids, for example). This method will work for arbitrary convex
+        geometries, but may result in more than one face for each truncated vertex.
+
+    Args:
+        poly (ConvexPolyhedron): Shape to truncate.
+        t (float):
+            Truncation depth as a fraction of initial edge length. Must be in [0.0, 0.5]
+        degrees (None | np.array | list, optional):
+            Vertex degrees to truncate. Defaults to None, which truncates all vertices.
+
+    Returns
+    -------
+        ConvexPolyhedron: Truncated polyhedron
+    """
+    new_vertices = _truncate(poly, t, degrees=degrees)
+    new_poly = ConvexPolyhedron(new_vertices)
+    new_poly.volume = poly.volume
+    return new_poly
+
+
+def dual(poly: ConvexPolyhedron):
+    """Dual polyhedron as defined by Conway, computed using halfspaces.
+
+    Args:
+        poly (ConvexPolyhedron): Shape to compute dual of
+
+    Returns
+    -------
+        ConvexPolyhedron: Dual polyhedron
+    """
+    dual = HalfspaceIntersection(
+        halfspaces=poly._simplex_equations,
+        interior_point=poly.centroid,
+        qhull_options="H Qt",
+    )
+    try:
+        dual_vertices = dual.dual_points[dual.dual_vertices]
+    except ValueError:
+        # If the halfspace intersection has double points, try manually
+        dual_vertices = np.unique(dual.dual_points, axis=1)
+    new_poly = ConvexPolyhedron(dual_vertices)
+    new_poly.volume = poly.volume
+    return new_poly
+
+
+def kis(poly: ConvexPolyhedron, k: float, degrees=ArrayLike):
+    """Pyramid-augment apolyhedron as defined by Conway's kis operator.
+
+    This method is implemented as the sequence of operators `dtd` to ensure a pyramid of
+    the proper height is raised on each face.
+
+    Args:
+        poly (ConvexPolyhedron):
+            Shape to pyramid augment.
+        k (float):
+            Pyramid height as a fraction of initial edge length. Should be in [0, 0.5].
+        degrees (None | ArrayLike, optional):
+            Face degrees to augment. Defaults to None.
+        normalize (bool, optional):
+            Whether to normalize output volume to match input. Defaults to True.
+
+    Returns
+    -------
+        ConvexPolyhedron: Augmented polyhedron polyhedron
+    """
+    new_poly = dual(ConvexPolyhedron(vertex_truncate(dual(poly), k, degrees=degrees)))
+    new_poly.volume = poly.volume
+    return new_poly

tests/test_shape_families.py

@@ -9,6 +9,7 @@
 from coxeter.families import (
     DOI_SHAPE_REPOSITORIES,
     ArchimedeanFamily,
+    CanonicalTrapezohedronFamily,
     CatalanFamily,
     Family323Plus,
     Family423,
@@ -29,6 +30,7 @@
 
 ATOL = 1e-15
 MIN_REALISTIC_PRECISION = 2e-6
+MIN_DECIMALS = 6
 MAX_N_POLY = 102
 TEST_EXAMPLES = 32
 
@@ -298,3 +300,18 @@ def test_new_pyramid_dipyramid():
             np.testing.assert_allclose(
                 shape.edge_lengths, comparative_shape.edge_lengths, atol=ATOL
             )
+
+
+@given(
+    integers(3, MAX_N_POLY),
+)
+def test_trapezohedra(n):
+    vertices = CanonicalTrapezohedronFamily.make_vertices(n)
+    poly = CanonicalTrapezohedronFamily.get_shape(n)
+    np.testing.assert_array_equal(vertices, poly.vertices)
+
+    assert np.isclose(poly.volume, 1.0)
+    assert len(np.unique(poly.edge_lengths.round(MIN_DECIMALS))) == 2 or np.allclose(
+        poly.edge_lengths, 1.0
+    )
+    assert all([len(face) == 4 for face in poly.faces])  # All faces are kites

tests/test_shapemoves.py

@@ -0,0 +1,83 @@
+# Copyright (c) 2015-2025 The Regents of the University of Michigan.
+# This file is from the coxeter project, released under the BSD 3-Clause License.
+
+import numpy as np
+import pytest
+
+from coxeter.families import ArchimedeanFamily, PlatonicFamily
+from coxeter.shapemoves import vertex_truncate
+
+
+def shapeeq(
+    poly1,
+    poly2,
+    test_vertices=False,
+    test_inertia_tensor=True,
+    test_surface_volume=True,
+    test_curvature=True,
+):
+    vx1, vx2 = poly1.vertices, poly2.vertices
+    assert poly1.num_vertices == poly2.num_vertices, (
+        f"Polyhedra do not have the same number of vertices "
+        f"({poly1.num_vertices} vs. {poly2.num_vertices})"
+    )
+
+    ex1, ex2 = poly1.edges, poly2.edges
+    assert np.shape(ex1) == np.shape(ex2), (
+        f"Polyhedra do not have the same number of edges "
+        f"({np.shape(ex1)[0]} vs. {np.shape(ex2)[0]})"
+    )
+
+    assert poly1.num_faces == poly2.num_faces, (
+        f"Polyhedra do not have the same number of faces "
+        f"({poly1.num_faces} vs. {poly2.num_faces})"
+    )
+
+    if test_inertia_tensor:
+        poly1.volume, poly2.volume = 1, 1
+        assert np.allclose(poly1.inertia_tensor, poly2.inertia_tensor), (
+            f"Inertia tensors do not match: {poly1.inertia_tensor.round(10)} vs. "
+            f"{poly2.inertia_tensor.round(10)}."
+        )
+
+    if test_vertices:
+        for vertex in vx1:
+            assert np.any([np.isclose(vertex, vertex2) for vertex2 in vx2]), (
+                f"Vertex {vertex} not found in poly2."
+            )
+
+    for edge in ex1:
+        assert edge in ex2 or edge[::-1] in ex2, f"Edge {edge} not found in poly2."
+
+    if test_surface_volume:
+        assert np.isclose(
+            poly1.volume / poly1.surface_area, poly2.volume / poly2.surface_area
+        ), (
+            f"{poly1.volume / poly1.surface_area} vs. "
+            f"{poly2.volume / poly2.surface_area}"
+        )
+
+    if test_curvature:
+        assert np.isclose(poly1.mean_curvature, poly2.mean_curvature)
+
+    return True
+
+
+@pytest.mark.parametrize(
+    "poly_name", ["Tetrahedron", "Cube", "Octahedron", "Dodecahedron", "Icosahedron"]
+)
+def test_truncation(poly_name):
+    uniform_truncation_depths = {
+        "Tetrahedron": 1 / 3,
+        "Cube": 1 / (2 + np.sqrt(2)),
+        "Octahedron": 1 / 3,
+        "Dodecahedron": 1 / (2 + (1 + np.sqrt(5)) / 2),
+        "Icosahedron": 1 / 3,
+    }
+    untruncated_poly = PlatonicFamily.get_shape(poly_name)
+    truncated_poly = ArchimedeanFamily.get_shape(f"Truncated {poly_name}")
+    test_poly = vertex_truncate(
+        untruncated_poly,
+        t=uniform_truncation_depths[poly_name],
+    )
+    assert shapeeq(truncated_poly, test_poly)