Welcome to antitile’s documentation!¶
Antitile is a package for manipulation of polyhedra and tilings using Python. It is designed to work with [Antiprism] but can be used on its own.
The bulk of this package pertains to Goldberg-Coxeter operations, as described in its own section.
Goldberg-Coxeter Operations on Polyhedra and Tilings¶
The Goldberg-Coxeter (GC) operation can be used to subdivide the faces of a polyhedron or tiling with triangular or square faces, by replacing the faces with a portion of a lattice of triangular or square faces. Geodesic subdivision is another name for the method on triangular faces; this method can be attributed to [Goldberg], [Coxeter], Fuller, or [Caspar] and Klug, depending on whether you’re a mathematician or a scientist :) This method often produces polyhedra with nice geometric qualities, for instance, local symmetry preservation, minimal distortion, etc.
Strictly, the GC operation defined in e.g. [Deza2004] and [Deza2015] is defined on a polyhedral graph, not a polyhedron. Embedding the graph into R3 produces a polyhedron. In the embedding, it is necessary to consider geometric questions like the exact placement of vertices that aren’t an issue with graphs; those are addressed here.
Notation¶
This text uses \(\Delta(a,b)\) for the triangular GC operator and \(\Box(a,b)\) for the quadrilateral GC operator. This contrasts with [Deza2004], who use \(GC_{k,l}(G_0)\) for both.
Contents¶
Master polygons¶
Lattices¶
The square and triangular lattices in the plane have the nice quality that all of their faces remain similar to each other with any affine transformation of the plane. This makes them a good candidate to use if, in a vague sense, we want to subdivide a face of a tiling or polyhedron into smaller faces that are somewhat similar to the original face.
The hexagonal lattice is dual to the triangular lattice. The square lattice’s dual is a lattice staggered with respect to the original lattice.
Relation to Gaussian and Eisenstein integers¶
The vertices of square and triangular lattices match with the elements of certain integral domains [1] in the complex plane. The Gaussian integers form a square lattice, and the Eisenstein integers form a triangular lattice. To make the parameterization of the \(\Delta(a,b)\) operator match the traditional parameterization used by [Goldberg] and [Coxeter], the Eisenstein integers can be paramaterized with \(u = e^{\pi i/3}\) instead of the usual \(w = e^{2\pi i/3}\). The relationship between the two is simple: \(a + b w = a - b + bu\).
| Type | Gaussian | Eisenstein |
|---|---|---|
| Form | \(a + b i\) | \(a + b u\) |
| Adjoined element | \(i =\sqrt{-1}\) | \(u = \frac{1}{2}(1 + i\sqrt 3) = e^{\pi i/3}\) |
| Units \((a, b)\) | \((1, 0) = 1\), \((0, 1) = i\), \((-1, 0) = -1\), \((0, -1) = -i\) | \((1, 0) = 1\) \((0, 1) = u\), \((-1, 1) = u^2 = u-1 = w\), \((-1, 0) = -1\), \((0, -1) = -u\), \((1, -1) = -u^2=1-u = -w\) |
| Algebraic Norm (Euclidean Function) \(|x|^2=x\overline x\) | \(a^2 + b^2\) | \(a^2 + ab + b^2\) |
| Multiplication \((a+bz) (c+dz)\) | \((ac-bd) + (bc+ad)i\) | \((ac-bd)+(bc+ad+bd)u\) |
| Division \((a+bz)/(c+dz)\) | \(\left({ac + bd \over c^2 + d^2}\right) + \left( {bc - ad \over c^2 + d^2} \right)i\) | \(\left({ac+ad+bd \over c^2 + cd+ d^2}\right) + \left({bc-ad \over c^2 + cd+ d^2}\right)u\) |
| Divisor t for prime test | 4 | 3 |
| Class sine | \(C = \sin^2(2\theta) = \frac{4a^2b^2}{(a^2+b^2)^2} = -\frac{(x^2-x^2)^2}{4|x|^4}\) | \(C = \sin^2(3\theta) = \frac{27a^2b^2(a+b)^2}{4(a^2+ab+b^2)^3} = -\frac{(x^3-x^3)^2}{4|x|^6}\) |
| Master polygon vertices | 0, \(x\), \(x(1+i)\), \(xi\) | 0, \(x\), \(xu\) |
| Master polygon center | \(x\frac{1+i}{2}\) | \(x\frac{1+u}{3}\) |
| Example master polygon (3,2) construction (center marked with white circle) |  |
 |
Elements \(x\) and \(y\) of a domain such that \(x = zy\), where \(z\) is a unit, are called associates. “Normal elements” are defined so that we have a “normal” form for associates: those elements having \(a > 0\) and \(b \ge 0\). All non-zero elements of the domain can be expressed as a normal element times a unit. All the units have \(|z| = 1\), so the associate elements can be thought of as the elements rotated around the origin (by steps of 90 degrees for the Gaussians, 60 degrees for the Eisenstein.)
An equivalence class can be defined from the associate relationship (excluding the zero element): all associated elements belong to the same equivalence class. The “normal element” is the preferred label \(x\) for that equivalence class.
Both these domains are Euclidean domains. That means that all non-zero elements have a unique factorization into prime elements and units, and what’s more, they can be factored using an extension of the Euclidean division algorithm on the natural integers. A prime element can only be divided by itself, an associate of itself, or a unit. More specifically, non-zero elements of these domains can be factored as:
where \(z\) is a unit, \(n_i\) is a natural integer, and \(p_i\) is a normal prime element of the domain.
An element of these domains is prime if:
- \(x = pz\), where \(z\) is a unit, \(p\) is a natural prime, and \(p = -1\mod t\), or
- \(|x|^2\) is a natural prime.
Relation between operators and master polygons¶
We define \(\Delta(a,b)\) for the triangular GC operator and \(\Box(a,b)\) for the quadrilateral GC operator. The GC operation can be thought of as taking each face of the polyhedron and replacing it with the corresponding master polygon described above, using the Gaussian integers for \(\Box(a,b)\) and Steineisen integers for \(\Delta(a,b)\). (The section on arbitrary surfaces will expand on this definition.). In fact, it has been shown that \(\Delta(a,b)\) corresponds to the equivalence class of associates in the Eisenstein integers, and \(\Box(a,b)\) corresponds to the equivalence class of associates in the Gaussian integers. There are some very nice consequences of this:
- Composition of operators corresponds to multiplication of complex numbers. If \(\Delta(a,b)\Delta(c,d) = \Delta(e,f)\), then \((a + bu)(c + du) = z(e + fu)\) for some unit z. A similar relation holds for \(\Box\).
- Since elements of these domains can be factored, operators can be “factored” into a sequence of operators. When \(x = a + bi\) is an element of the Gaussian integers, and \(x = z p_1^{n_1}p_2^{n_2} \cdots p_k^{n_k}\), then \(\Box(a,b) = \Box^{n_1}(p_1)\Box^{n_2}(p_2)\cdots\Box^{n_k}(p_k)\), and similarly for the Eisenstein integers and \(\Delta(a,b)\). It also makes sense to talk about “prime” operators that can’t be decomposed further.
\(\Delta(1,0)\) and \(\Box(1,0)\) are identity operators: where defined, they cause no change to the subdivision.
In the terminology of geodesic domes,
- Class I operators have \(b=0\).
- Class II operators have \(b=a\).
- Class III operators are all others, and occur in chiral pairs. \(\Delta(a,b)\) is the chiral pair of \(\Delta(b,a)\), and the same for \(\Box(a,b)\) and \(\Box(b,a)\).
Here we also introduce the ‘’class sine’‘, to have a smooth measure of how much like a class I or class II operator a given operator is. If an operator is Class I, \(C=0\); if Class II, \(C=1\), and Class III is in between.
Chiral pairs correspond to conjugation of the complex number \(a+bi\) or \(a+bu\). That is, there exists some unit z such that \(z(a-bi) = b+ai\), or \(z(a+b\bar{u}) = b+au\). Interestingly, the notation used in Conway notation for the chiral pair of an operator is an overbar, also used to denote complex conjugation. A consequence of this relation between chiral pairs is that the composed operator \(\Delta(a,b)\Delta(b,a)\) will always be a Class I operator, since a complex number times its conjugate is a real number. For every Class III \((a,b)\), there will exist three operators \(\Delta(a,b)\Delta(b,a)\), \(\Delta(a,b)\Delta(a,b)\), and \(\Delta(b,a)\Delta(b,a)\) with the same algebraic norm. The same holds for \(\Box(a,b)\).
Footnotes
| [1] | A type of ring. |
Goldberg-Coxeter Operation on arbitrary tilings and polyhedra¶
This section will discuss tilings on arbitrary surfaces before we discuss the
spherical case, i.e. the behavior of the program with options
gcopoly -n -p=flat. Usually the GC operation is described on a usual
polyhedra of genus 0, but it can be applied to most surfaces. The main caveat
is that class III operators doesn’t work on non-orientable surfaces: because
the operators are chiral, they depend on an underlying orientable surface.
Definition of GC operation¶
The GC operation comprises these steps:
- Choose a master polygon corresponding to \(\Box(a,b)\) or \(\Delta(a,b)\).
- Replace each face of the base polyhedron with a stretched, transformed version of the master polygon that covers the base face.
- Stitch together any edges that might cross the edges of the base face.
There are some places in the literature where the GC operation is defined as the dual of this operation. Particularly, in the fullerene literature, most writers are interested in hexagonal faces, so they use the dual.
Coordinate form¶
To transform the master polygon to the target face polygon, this program generally uses a standardized coordinate form.
Triangles can be are specified by barycentric coordinates \(\beta_i\) where \(\beta_1 + \beta_2 + \beta_3 = 1\), and then the corresponding vertex is given by \(\mathbf v = \beta_1\mathbf v_1+\beta_2\mathbf v_2+\beta_3\mathbf v_3\). Given \(\mathbf v\) and \(\mathbf v_i\), \(\beta_i\) can be found by using e.g. the Moore–Penrose inverse to solve the system given by \(\beta_1 + \beta_2 + \beta_3 = 1\) and \(\mathbf v = \beta_1\mathbf v_1+\beta_2\mathbf v_2+\beta_3\mathbf v_3\).
Quadrilaterals are instead specified by what we’ll call “xy coordinates” where \(x\) and \(y\) are between 0 and 1. That transformation is:
Unlike triangles, quadrilaterals may have points that do not share a common plane: they may be skew quadrilaterals. If the quadrilateral is a skew quadrilateral, x and y smoothly parameterize a surface over that skew quadrilateral. In this program, we don’t need to find the xy coordinates for a point in an arbitrary quadrilateral. There are some occasions where the xy coordinates are determined for a square in the plane, in which case they can be converted by rotation and scaling.
Linear index form¶
Another way to think about the master polygon is to count lines in the plane that cross it. This method can break down at the base vertices, but we already know what those are supposed to be, so that’s not a practical problem.
For \(\Delta(a,b)\), draw an equilateral triangle in the plane. Along one edge, mark a+b+1 evenly spaced points (including the base vertices). Going clockwise, mark a+1 points along the next edge, then b+1 along the last. Number the points from 0 to a+b along the first side, then along the remaining two sides sequentially (counting the point shared between the other two sides only once). Draw lines between points with the same number. Then take those lines and make copies rotated 120° and 240° about the center. The vertices of the master triangle correspond to the points where 3 lines intersect (not necessarily the marked points), and the linear index is a 3-tuple of the number we gave to the original line, the 120° line, and the 240° line that intersect at that vertex.
For \(\Box(a,b)\), draw a square in the plane. Mark b+1 vertices along the left edge of the square (including the base vertices) and a+1 vertices along the bottom. Number the points from 0 to a+b, starting at the upper left. Copy these points, rotate them by 180° about the center, and reverse their numbering. Draw a line between each point with the same number. Then take those lines and make a copy rotated by 90° about the center. The vertices of the master square correspond to the points where 2 lines intersect (not necessarily the marked points), and the linear index is a 2-tuple of the number we gave to the original line and the rotated line that intersect that that vertex.
Given the Gaussian or Eisenstein coordinates of the vertices in the master polygon, they can be directly converted to linear indexes.
Triangular: When the vertex is \(c + du\), in normal form, and the master triangle is given by (a, b), then the linear indexes satisfy \(\sum \ell_i = 2a+b\), with
Quadrilateral: When the vertex is \(c + di\), in normal form, and the master square is given by (a, b), then
The primary use for linear index in the gcopoly program is “stitching”
neighboring polygons together at their edges. This operation is easy to
do by eye but a little complicated to automate. Linear indexes have the
advantage that they are integers, not floats, and therefore rounding error
isn’t an issue and the stitching can be done exactly.
This program, when creating the master polygon, leaves in some vertices that are outside the polygon but share edges with vertices inside it. This is done in part to make stitching easier: the program determines which vertices in adjacent faces match up, and collapses them into one. In detail, the stitching procedure is as follows:
- Select two faces that share an edge. Make a copy of the linear indexes for each.
- Rotate the order of the linear index of both faces so that the vertices shared by the edge have the same lindex. (If a class I or II operation, can reverse the tuple order too to accomodate unoriented surfaces.)
- On one face, replace the linear index with \(\mathbf o - \mathbf \ell\), where \(\mathbf o = (a+b, a, 2a+b)\) for \(\Delta(a,b)\) and \(\mathbf o = (a+2b, b)\) for \(\Box(a,b)\).
- Collapse vertices on different faces with the same linear index.
Surfaces with boundaries present a problem because there are tiles that don’t have a neighboring tile for each edge. In the language of [Deza2004], the skeleton graph of this tiling has vertices with degree 2; that paper doesn’t even define their operation on such graphs. As mentioned earlier, this program leaves in vertices that are outside the master polygon but are connected to vertexes inside it. Correspondingly, the program will retain those vertices when it transforms the master polygon. This means that surfaces with boundaries may have faces outside the original boundary if a Class II or III operator is applied.
(However, see the appendix gco-a-improper for information on the improper spherical tilings that may have graphs with vertices of degree 2 or less. These are not surfaces with boundaries, and can be handled differently.)
Composition: Topology vs Geometry¶
On the plane, equality of operators holds in both a topological and geometric sense: if \(\Box(a,b)\Box(c,d) = \Box(e,f)\), then the left-hand side results in the same placement of vertices on the plane. This topological equality carries over to tilings on arbitrary surfaces, but not necessarily the geometric sense. Furthermore, commutativity and associativity hold for topology, but not necessarily geometry. So considering the exact position of vertices in space, where \((a,b)\), \((c,d)\), and \((e,f)\) are complex numbers in the appropriate space:
- \(\Delta(a,b)\Delta(c,d)\) may not equal \(\Delta(e,f)\), where \((e,f)\) is the normal form of the product of \((a,b)\) and \((c,d)\)
- \(\Delta(a,b)\Delta(c,d)\) may not equal \(\Delta(c,d)\Delta(a,b)\)
- \(\Delta(a,b)(\Delta(c,d)\Delta(e,f))\) may not equal \((\Delta(a,b)\Delta(c,d))\Delta(e,f)\)
and the same for \(\Box\).
This fact actually turns out to be useful, because it gives us some flexibility in vertex placement.
Goldberg-Coxeter Operation on spherical polyhedra¶
Spherical (genus 0) polyhedra are the standard use case for the GC operation. There are many methods in the geodesic dome literature for placing the vertices on the surface of a sphere. Many of these are synthetic (step-by-step geometric contstructions) rather than analytic (expressed in equations). For a computer application, analytic forms can be more convenient. This section might not cover every method ever written, but it expresses the major methods in analytic form and introduces a couple new ones.
In this section, we’ll use \(\mathbf v^*\) to denote a pre-normalized vector:
Gnomonic¶
gcopoly -p=flat
The gnomonic projection was known to the ancient Greeks, and is the simplest of the transformations listed here. It has the nice property that all lines in Euclidean space are transformed into great circles on the sphere: that is, geodesics stay geodesics. This is in fact the motivation for the name “geodesic dome”: Buckminster Fuller used this projection to project triangles on the sphere. This is referred to as Method 1 in geodesic dome circles. The main downside is that the transformation causes shapes near the corners to appear bunched up; this is particularly bad for larger faces e.g. on the tetrahedron.
In general, the gnomonic projection is defined as:
- To sphere: \(\hat{\mathbf v} = \frac{\mathbf p}{\|\mathbf p\|}\)
- From sphere: \(\mathbf p = \frac{r\hat{\mathbf v}} {\hat{\mathbf n} \cdot \hat{\mathbf v}}\)
where \(\mathbf p\) is a point on a plane given in Hessian normal form by \(\hat{\mathbf n} \cdot \mathbf p = r\). Projection from Euclidean space to the sphere is literally just normalizing the vector.
For the Goldberg-Coxeter operation, this amounts to just normalizing the vectors produced by the coordinate form:
where \(\beta_i\) are (planar) barycentric coordinates and \(x,y\) are x-y quadrilateral coordinates.
The case of barycentric coordinates can be thought of as generalized barycentric coordinates, where \(\mathbf v = \sum\beta^\prime_i\mathbf v_i\) still holds but \(\sum \beta^\prime_i\) is not necessarily =1. If the generalized coordinates are \(\beta^\prime_i\), then \(\beta^\prime_i = \frac{\beta_1} {\beta_1 \mathbf v_1 + \beta_2 \mathbf v_2 + \beta_3 \mathbf v_3}\). On the interior of the triangle, \(\sum \beta^\prime_i > 1\).
Spherical areal¶
gcopoly -p=other
This method applies to triangles only. The above section discussed generalized barycentric coordinates, which removes the restriction that \(\sum \beta_i = 1\). Instead we look to the relation between barycentric coordinates and area; \(\beta_i\) is the proportion of area in the triangle that is opposite the vertex \(v_i\). Let \(\Omega\) be the spherical area (solid angle) of the spherical triangle and \(\Omega_i = \beta_i\Omega\) be the area of the smaller triangle opposite vertex \(v_i\).
While this method introduces less bunching-up than the gnomonic method, the equation to find \(\hat{\mathbf v}\) given \(\beta_i\) is much more complicated than for barycentric coordinates. Let \(h_i = \sin\Omega_i\left(1+\mathbf v_{i-1}\cdot\mathbf v_{i+1}\right)\), and \(\mathbf g_{i} = \left(1+\cos \Omega_{i}\right) \mathbf v_{i-1} \times \mathbf v_{i+1} - \sin\Omega_{i}\left(\mathbf v_{i-1} + \mathbf v_{i+1}\right)\) where the subscripts loop around: 0 should be interpreted as 3 and 4 should be interpreted as 1. Then
such that \(\mathbf G \hat{\mathbf v} = \mathbf h\). To clarify, \(\mathbf G\) is the 3x3 matrix where the ith column is \(\mathbf g_i\), and \(\mathbf h\) is the column vector where the ith element is \(h_i\). The vector \(\hat{\mathbf v}\) can be solved for using standard matrix methods. (This mess can be derived from the formula for solid angle in the appendix.)
Naive Slerp¶
gcopoly -p=nslerp, gcopoly -p=nslerp2
This method shares with the gnomonic method an analytic form for the transformation from Euclidean space to the sphere, but doesn’t have the problem with bunching up. It has two downsides: there is no analytic form for the reverse transformation, and it can only be used on equilateral faces.
The Naive Slerp method on a triangular face resembles a naive extension of spherical linear interpolation (Slerp) to barycentric coordinates, thus the name. The Naive Slerp methods reduce to slerp on the edges of the face.
Let \(\cos(w) = \mathbf v_i \cdot \mathbf v_{i+1}\) for all \(i\). (As usual, the subscripts loop around.)
Triangle:
Quadrilateral 1:
- \(\gamma_1 = (1-x)(1-y)\),
- \(\gamma_2 = x(1-y)\),
- \(\gamma_3 = xy\),
- \(\gamma_4 = (1-x)y\)
Quadrilateral 2:
- \(s_1 = \sin (w(1-x))\sin (w(1-y))\),
- \(s_2 = \sin (wx)\sin (w(1-y))\),
- \(s_3 = \sin (wx)\sin (wy)\),
- \(s_4 = \sin (w(1-x))\sin (wy)\)
Because the projected edges already lie on the sphere, there is a lot of freedom in how to adjust \(\mathbf v^*\) to lie on the sphere. The easiest is just to centrally project the vertices, that is, to normalize \(\mathbf v^*\) like we have been. Another option is to perform a parallel projection along the face normal. (See appendix for a formula for the “normal” to a skew face.) We need the parallel distance p from the vertex to the sphere surface in the direction of the face normal \(\hat{\mathbf n}\), such that \(\hat{\mathbf v} = \mathbf v^* + p\hat{\mathbf n}\). p is given by:
p can also be approximated as \(\widetilde{p} = 1 - \|\mathbf v^*\| \leq p\), which is somewhat fewer operations and doesn’t require calculation of the face normal.
Sometimes the best polyhedron comes from a compromise of the central and parallel projections. Choose a constant k, typically between 0 and 1, then:
p may be replaced by \(\widetilde{p}\).
If our goal is to optimize a measurement of the polyhedron, we can do a 1-variable optimization on k, which is usually more tractable than the multivariate optimization of the location of every vertex.
(Technically, you can project in almost any direction, not just that of the face normal, but most other choices don’t produce anything remotely symmetric.)
Great circle intersection¶
gcopoly -p=gc
This method draws great circles between points on the edges of the polygon, and uses the points of intersections of those great circles to determine the vertices. This is the only method that directly uses the linear indexes defined in the previous chapter. In the geodesic dome world, this is called Method 2, although this description is considerably more complicated to accomodate Class III grids and quad faces.
Specify the linear index as \(\ell = (e,f,g)\) for triangular faces or \(\ell = (e,f)\) for quad faces. Using slerp, calculate the points \(\mathbf{\hat{b}}_{i,j,k}\) where each line from the breakdown structure crosses the face edge. i is the coordinate of the linear index, j is the linear index (\(\ell_i = j\)), and k is which point of intersection with the polygon. The line corresponds to a great circle normal given by \(\mathbf{\hat{n}}_{i,j} = \frac{\mathbf{\hat{b}}_{i,j,0} \times \mathbf{\hat{b}}_{i,j,1}}{\|\dots\|}\).
We are going to calculate the intersection of these great circle normals. The intersection of two planes is a line: the intersection of two great circles is two antipodal points. We need to choose the point on the correct side of the sphere. Let \(\mathbf{c}\) be the centroid of the face: then \(\mathbf{v}\) is on the right side of the sphere if \(\mathbf{v} \cdot \mathbf{c} >0\). If not, just multiply \(\mathbf{v}\) times -1 to put it on the right side.
For quad faces, there are only two intersecting great circles, so the new vertices are \(\mathbf{v^*}_{\ell} = \mathbf{\hat{n}}_{1,\ell_{1}} \times \mathbf{\hat{n}}_{2,\ell_{2}}\) (possibly times -1).
For triangular faces, there are three intersecting great circles, and unlike on the plane, on the sphere they do not necessarily intersect in the same place. Each pair of great circles forms a vertex of a triangle as \(\mathbf{\hat{v}}_{\ell, m} = \frac{\mathbf{\hat{n}}_{m,\ell_{m}} \times \mathbf{\hat{n}}_{m+1,\ell_{m+1}}}{\|\dots\|}\), Make sure all the \(\mathbf{\hat{v}}_{\ell, m}\) lie on the correct side of the sphere, and then take the centroid of that triangle to get the vertex: \(\mathbf{v^*}_{\ell} = \sum_m \mathbf{v}_{\ell, m}\). (It is not strictly necessary to take the centroid: the sum of the unnormalized \(\mathbf{v}_{\ell, m}\) will also be a point somewhere within the triangle.)
Summary of methods¶
| Method | Gnomonic | Spherical areal | Naive slerp | Great circle intersection |
|---|---|---|---|---|
| Geodesic dome name | Method 1 | New | New | Method 2 |
| Input | Coordinates (barycentric or xy) | Barycentric coordinates | Coordinates (barycentric or xy) | Linear index (triangular or quadrilateral) |
| Adjustment to sphere | Central projection | Not needed | Any projection | Central projection |
| Face | \(\Omega < 2\pi\) | \(\Omega < 2\pi\) | \(\Omega <= 2\pi\), must be equilateral | \(\Omega < 2\pi\) |
Multi-step methods¶
As mentioned earlier, the operators \(\Delta(a,b)\) and \(\Box(a,b)\) may be able to be decomposed into a series of smaller operators. Many of the smaller operators are constrained by symmetry: in particular, \(\Delta(2,0)\) adds vertices at the midpoints of each edge, independent of the method used. Method 3 in geodesic dome terminology is simply repeated application of \(\Delta(2,0)\).
In a more general sense, an operator can be factored into a series of “prime” operators, and applied in order from small to large. The faces of the polyhedron will become progressively smaller and therefore progressively flatter, and as the faces get flatter, the differences between methods becomes smaller. As an example, \(\Box(16,4) = \Box^4(1,1)\Box(4,1)\), so apply the highly-symmetric operator \(\Box(1,1)\) (which creates one vertex at the centroid of a face) four times and then \(\Box(4,1)\) once with a simple method like Gnomonic.
geodesic in Antitile performs class II and III subdivision by finding the
smallest class I operator that can be decomposed into the desired operator
and some other factor. Effectively, given \(\Delta(a,b)\), it finds the
smallest n such that \(\Delta(a,b)\Delta(c,d) = \Delta(n,0)\) for some
c and d. That is given by \(n = \frac{t(a,b)}{\gcd(a,b)}\),
where gcd is the greatest common divisor.
It calculates \(\Delta(n,0)\) using method 2 and then uses the
vertices from that that are shared with \(\Delta(a,b)\).
Skew faces¶
Skew faces are impossible on a polyhedra with triangular faces. On a polyhedron with quadrilateral faces, however, all of the above methods produce skew faces. There are basically two solutions to the issue. The first is to treat the polyhedron purely as a spherical polyhedron: all the faces are curved tiles on the surface of a sphere, and we can ignore whether they’re skewed in Euclidean space. The second is to canonicalize the polyhedron. As per [Hart1997], all convex polyhedra can be put into a unique canonical form such that:
- All the edges are tangent to the unit sphere,
- The origin is the average of the points at which the edges touch the sphere, and
- The faces are flat (not skew)
The canonical program in Antiprism performs canonicalization via a simple
iterative process. The vertices of the faces probably do not lie on the
unit sphere. If a polyhedron created by Goldberg-Coxeter
operations is to be canonicalized, the choice of method does not matter except
as a starting point.
Choosing a method¶
Which method is better depends on your criteria. If your only criteria is speed, gnomonic is the simplest and fastest to run, and produces reasonable results on typical cases like the icosahedron and cube.
A common criteria is to have a polyhedron that minimizes some measure of deviation. That could be a number of things:
- Thomson energy
- Edge length (spherical or euclidean)
- Aspect ratio for each face (spherical or euclidean)
- Face area (spherical or euclidean, only well-defined for non-skew faces)
- Face bentness (only applies to quad faces)
Often the objective will be to minimize either the maximum value or the ratio
of maximum to minimum values. This package includes, in the data directory,
a csv file containing measurements of polyhedra produced by GC operations
on the icosahedron, octahedron, tetrahedron, and cube. No method always comes
out the winner, but some patterns emerge:
- Gnomonic doesn’t ourperform any other method (except on speed)
- Naive slerp often outperforms all other methods. This is true when the k-factor is fixed at 1 and moreso when k-factor is optimized.
- For Class I grids, all methods except Gnomonic and Areal produce optimal edge lengths on reasonably small faces.
- Great circle works very well for Class I grids, especially with the tweak mentioned above.
- On quad faces, the first Naive Slerp generally outperforms the second.
geodesicin Antiprism performs well with regards to edge length ratio.- Spherical areal performs well for certain measures on large faces, e.g. edge length ratio and face area for the tetrahedron
- No method (aside from canonicalization) consistently creates all perfectly planar quadrilateral faces.
Due to symmetry, all methods produce the same results for parameters (1,0), (1,1), and (2,0). With \(\Delta(3,0)\), great circle and naive slerp (with any k value) produce the same results.
[Altschuler] conjectures that for a geodesic sphere is to being class I, the lower its Thomson energy will be. This seems to be true for most other measurements, and quadrilateral faces: given GC operators with the same norm, the operators with the lower class sine seem to produce lower Thomson energy, edge length ratios, differences in aspect ratio, face area ratios, and so forth.
To demonstrate, here is a table with metrics for \(\Delta(5,3)I\) and
\(\Delta(7,0)I\) for a number of methods. Both have a T-value of 49, 492
vertices, and 980 faces. k values other than 1 are chosen to optimize one of
the metrics. E is the Thomson energy of the vertices: Eo = 115005.2559 is
approximately the optimal Thomson energy for 492 vertices, found using
repel in Antiprism. “Edge ratio” is the length of longest edge divided by
the smallest edge, minus one. “Face ratio” is the area of the largest face
divided by the smallest face, minus one.
| a | b | Method | k | E - Eo | Edge ratio | Face ratio |
|---|---|---|---|---|---|---|
| 7 | 0 | gnomonic | 134.91 | 36.42% | 69.66% | |
| 7 | 0 | areal | 11.00 | 17.57% | 20.26% | |
| 7 | 0 | antiprism geodesic (same as great circle) | 4.40 | 17.19% | 13.42% | |
| 7 | 0 | great circle (tweak) | 2.14 | 17.19% | 10.85% | |
| 7 | 0 | nslerp | 1.00 | 2.91 | 17.19% | 11.74% |
| 7 | 0 | nslerp | 1.04 | 2.79 | 17.19% | 11.49% |
| 7 | 0 | nslerp | 1.26 | 2.47 | 17.19% | 11.71% |
| 7 | 0 | nslerp (tweak) | 1.00 | 3.09 | 17.19% | 11.84% |
| 7 | 0 | nslerp (tweak) | 1.05 | 2.91 | 17.19% | 11.50% |
| 7 | 0 | nslerp (tweak) | 1.30 | 2.53 | 17.19% | 11.74% |
| 5 | 3 | gnomonic | 154.96 | 43.30% | 87.17% | |
| 5 | 3 | great circle | 29.45 | 25.07% | 38.81% | |
| 5 | 3 | great circle (tweak) | 27.30 | 25.11% | 36.61% | |
| 5 | 3 | areal | 17.06 | 22.69% | 30.34% | |
| 5 | 3 | antiprism geodesic | 8.14 | 21.44% | 22.74% | |
| 5 | 3 | nslerp | 0.49 | 11.01 | 20.45% | 25.90% |
| 5 | 3 | nslerp | 1.00 | 5.50 | 22.25% | 19.20% |
| 5 | 3 | nslerp | 1.50 | 3.75 | 24.32% | 18.25% |
| 5 | 3 | nslerp | 1.74 | 4.11 | 25.23% | 17.87% |
| 5 | 3 | nslerp (tweak) | 0.51 | 11.03 | 20.55% | 25.84% |
| 5 | 3 | nslerp (tweak) | 1.00 | 6.01 | 22.28% | 19.75% |
| 5 | 3 | nslerp (tweak) | 1.55 | 4.08 | 24.55% | 18.76% |
| 5 | 3 | nslerp (tweak) | 1.78 | 4.40 | 25.47% | 18.47% |
Convexity¶
As a side note, often the geometries produced by the GC operation on a triangle-faced seed are convex polyhedra. In that case, the stitching operation outlined in the last section can be replaced with the operation of taking the convex hull of vertices. Almost all GC operations on an icosahedron or octahedron under all methods and reasonable values for k produce a convex polyhedron.
Appendix: Spherical (and other) geometry with vectors¶
In the geometric dome world, geometry is usually synthetic rather than analytic. That is, geometric constructions are usually expressed in terms of step-by-step geometric constructions rather than in equations and numeric values. See, for instance, [Kenner]. Implementing a synthetic geometric construction on a computer can be a pain; analytic specifications are more amenable to programming idioms.
There are a number of ways to numerically specify points on a sphere. By far the most common is by latitude and longitude, which appear on every modern map of the Earth and probably some other planets. Latitude and longitude are familiar and convenient, but performing extended geometry calculations using latitude and longitude is a complicated task.
Doing spherical geometry using a 3-component unit vector is more convenient in a number of ways: the equations are often simpler, and there are no singularities at the poles. This appendix isn’t a full course in spherical geometry with vectors, but is here to explain the notation and methods that are used in this text.
A 3d unit sphere can be defined as the set of all unit vectors in 3-space; i.e., vectors \(\mathbf v = [v_x, v_y, v_z]\) such that the vector norm \(\|\mathbf v \|=1\). Unit vectors are often denoted using a hat: \(\hat{\mathbf v}\). We’ll often find ourselves normalizing vectors, so if the numerator and denominator are the same, we’ll supress the denominator with an ellipsis like so:
Basics¶
The shortest distance between two points in Euclidean space is a straight line. On the sphere, the shortest distance is an arc of the great circle between those points. The great circle is the intersection of the sphere and a plane passing through the origin. A plane through the origin can be specified as \(\hat{\mathbf n} \cdot \mathbf v = 0\), where \(\hat{\mathbf n}\) is a unit vector normal to the plane; this vector \(\hat{\mathbf n}\) can be used to specify a great circle. Given two points \(\mathbf{\hat{v}_1, \hat{v}_2}\) on the sphere, the \(\hat{\mathbf n}\) of the great circle between those two points is (up to normalization) their cross product: \(\mathbf{\hat{n}} = \frac{\mathbf{\hat{v}}_1 \times \mathbf{\hat{v}}_2}{\|\dots\|}\)
Small circles are the intersection of the sphere with a plane not through the origin. All planes may be specified in Hessian normal form as \(\mathbf{\hat{n}} \cdot \mathbf v = r\), where r is the minimum distance between the plane and the origin.
| r | Intersection with unit sphere |
|---|---|
| 0 | Great circle |
| > 0 and < 1 | Small circle |
| 1 | Point |
| > 1 | None |
Measurements¶
| Measurement | Euclidean | Spherical |
|---|---|---|
| Length | \(\eta = \|\mathbf v_1-\mathbf v_2\|\) | Central angle: \(\theta = \arctan\left( \frac{\|\mathbf{\hat{v}}_1 \times \mathbf{\hat{v}}_2\|} {\mathbf{\hat{v}}_1 \cdot \mathbf{\hat{v}}_2}\right)\) |
| Angle | \(\cos \phi_1 = \frac{(\mathbf v_1 - \mathbf v_2) \cdot (\mathbf v_1 - \mathbf v_3)} {\|\mathbf v_1 - \mathbf v_2\|\|\mathbf v_1 - \mathbf v_3\|}\) | Dihedral angle: \(\mathbf{\hat{c}}_{12} = \frac{\mathbf{\hat{v}}_1 \times \mathbf{\hat{v}}_2}{\|\dots\|}\), \(\mathbf{\hat{c}}_{13} = \frac{\mathbf{\hat{v}}_1 \times \mathbf{\hat{v}}_3}{\|\dots\|}\), \(\cos\phi_1 = \mathbf{\hat{c}}_{12} \cdot \mathbf{\hat{c}}_{13}\) |
| Triangle area | \(A = \frac{\|(\mathbf v_1-\mathbf v_3)\times (\mathbf v_2-\mathbf v_3)\|}{2}\) | Solid angle: \(\tan(\Omega/2) = \frac{|\mathbf{\hat{v}_1} \cdot \mathbf{\hat{v}}_2 \times \mathbf{\hat{v}}_3|} {1+\mathbf{\hat{v}}_1\cdot \mathbf{\hat{v}}_2+\mathbf{\hat{v}}_2 \cdot \mathbf{\hat{v}}_3+\mathbf{\hat{v}}_3\cdot \mathbf{\hat{v}}_1}\) [Oosterom] |
In general, the spherical measures approach the Euclidean measures when the measures are small. (A sphere is locally Euclidean.)
Constructions¶
| Construction | Euclidean | Spherical |
|---|---|---|
| Mean (midpoint when n=2, centroid when n=3) | \(\mathbf v_\mu = \frac{\sum\mathbf v_i}{n}\) | \(\mathbf{\hat{v}}_\mu = \frac{\sum\mathbf{\hat{v}}_i}{\|\dots\|}\) |
| Interpolation | \(\mathrm{Lerp}(\mathbf{v_1}, \mathbf{v_2}; t) = (1-t) \mathbf{v_1} + t \mathbf{v_2}\) | \(\mathrm{Slerp}(\mathbf{\hat{v}_1}, \mathbf{\hat{v}_2}; t) = \frac{\sin {((1-t)w)}}{\sin (w)} \mathbf{\hat{v}_1} + \frac{\sin (tw)}{\sin (w)} \mathbf{\hat{v}_2}\) where \(\cos(w) = \mathbf{\hat{v}_1} \cdot \mathbf{\hat{v}_2}\) |
Face normals¶
This program uses this definition for the normal to a Euclidean polygon:
i should be treated as if it’s \(\mod n\): it loops around. This definition allows for a somewhat sensible extension to skew polygons: the normal points in the general direction that’s expected.
The normal will be outward-facing if the points are ordered counterclockwise, and inward-facing if the points are ordered clockwise.
Face bentness¶
There’s no standard measure of face bentness, so this program uses an ad-hoc measure that seems to work well. This program measures the bentness of a face with 4 or more vertices by these steps:
- Let \(\mathbf x_i = \mathbf{v}_i - \bar{\mathbf{v}}\), where \(\bar{\mathbf{v}}\) is the (Euclidean) average of the points
- Calculate the SVD decomposition of the matrix that has \(\mathbf x_i\) as rows (or columns). We only need the singular values: since we’re in 3d space, there will be 3 singular values.
- The “bentness” is the smallest singular value divided by the sum of the other two singular values.
Appendix: Goldberg-Coxeter Operation on improper spherical polyhedra¶
There are some improper spherical polyhedra that cannot be tesselated with the traditional geodesic dome methodology but can with some methods developed earlier. The Naive Slerp methods are able to project a triangle or square onto a full hemisphere, including the boundary.
Transformations to the disk¶
When the base vertices all lie on the same great circle (i.e. in a plane through 0), Naive Slerp before projection or normalization transforms the triangle or square to the disk. Unlike Triangular and Quadrilateral #1 Naive Slerp, Quadrilateral #2 simplifies nicely on the disk. If the vertices of the quadrilateral are \(\left<\pm 1, 0, 0\right>\) and \(\left<0, \pm 1, 0\right>\), and \(x,y\) are the xy coordinates, then:
The radial coordinate r of this vector in polar form has a nice form. (The angle coordinate \(\theta = \operatorname {arctan2} (v_y,v_x)\) doesn’t simplify usefully.)
There are also some disk-specific transformations, accessible with
gcopoly -p=disk. These amount to converting the vertex to polar form
and scaling the r variable. Unlike Naive Slerp, there are points where the
transformation is not smooth.
Triangular:
Quadrilateral:
Dihedra¶
The 3-dihedron
The 4-dihedron
With the transformations defined above, performing the GC operation on a dihedron is fairly straightforward. These disk operations can be projected with a k factor the same as naive slerp, although here we require k to be >0: otherwise, all points would be projected onto the great circle.
The only caveat has to do with the presence of vertices with valence 2. On the quadrilateral, the valence-2 vertex and its adjacent edges can be smoothed into a single edge. On the triangle, this produces dangling faces, which can be removed. In both cases, one can also add or switch edges so that the vertex is no longer degenerate.
The process of doing the GC operation on a dihedron resembles stitching together two flat sheets and inflating them, like the children’s novelty the whoopi cushion. The dual of that could be called a Whoopi Goldberg polyhedron. (Thank you, I’ll be here all night.)
Monohedron¶
The monohedron.
The spherical monohedron is the “polyhedron” with one face, one vertex, and no edges: the entire sphere is one big face. (Hey, it satisfies the Euler equations.) Topologically, we can tile this surface by taking the disk defined earlier and collapsing the boundary to a single point. This amounts to projecting the disk in the manner of a map projection of the entire Earth, which is a well-studied problem. Given \(\theta = \arctan2(v_x, v_y)\) and \(\mathbf{\hat{v}} = (\sin(\phi) \cos(\theta), \sin(\phi) \sin(\theta), \cos(\phi))\), two useful projections used in cartograph can be described as such:
Lambert azimuthal equal-area:
Azimuthal equal-distance:
Really any continuous, monotonic function of \(\|\mathbf v\|\) that maps 0 to 0 and 1 to \(\pi\) could be used to calculate \(\phi\).
We’ll call these “balloon polyhedra” by analogy with blowing up a balloon.
balloon.py can be used to calculate these. Because of the radical change
to the boundary, many small-parameter GC
operations produce the same polyhedron. Quad faces may be reduced to triangle
faces, and triangle faces may be reduced to degenerate faces. Many produced
polyhedra are not convex. Class III polyhedra tend to have dangling faces.
Appendix: Towards Goldberg-Coxeter operation on mixed polyhedra¶
The methods discussed up to now only handle polyhedra with a single type of face, either 3- or 4-sided. This section discusses the extension of the Goldberg-Coxeter operation to polyhedra with mixed triangle and quad faces. It’s possible to divide a polyhedron by applying a certain \(\Delta\) operator to the triangle faces, a certain \(\Box\) operator to the quadrilateral faces, and dealing with what happens over the edges of the original polyhedron.
Two operators are ‘’compatible’’ if:
- The same number of vertices lie along the base edge.
- The same number of edges cross the base edge, excluding ones that meet a vertex at the edge.
- The same number of edges lie on the base edge.
- Edges that cross the base edge do not create digons (faces with 2 edges).
For the operators \(\Delta(a,b)\) and \(\Box(c,d)\), the first 3 requirements can be distilled to these equations:
- \(\gcd(a, b) = \gcd(c, d) = g\)
- \(c + d + g = 2(a + b)\). This equation has no solution if c and d are both odd numbers.
The 4th requirement is more complicated to put in equational form.
Consequences of these rules include:
- \(\Box(1,0)\) is compatible with \(\Delta(1,0)\)
- Iff \(\Box(a,b)\) is compatible with \(\Delta(c,d)\), then \(\Box(na,nb)\) is compatible with \(\Delta(nc,nd)\) for all positive integers n.
- \(\Box(n,0)\) is compatible with \(\Delta(n,0)\) (Class I)
- Iff \(\Box(a,b)\) is compatible with \(\Delta(c,d)\), then \(\Box(b,a)\) is compatible with \(\Delta(d,c)\).
- \(\Delta(n,n)\) is compatible with \(\Box(n,2n)\) and \(\Box(2n,n)\) (Class II triangles)
- \(\Box(n,n)\) is not compatible with any \(\Delta\). (Class II squares)
- If \(\Box(a_1,b_1)\) is compatible with \(\Delta(c_1,d_1)\), and \(\Box(a_2,b_2)\) is compatible with \(\Delta(c_2,d_2)\), it does not follow that the composed subdivision \(\Box(a_1,b_1)\Box(a_2,b_2)\) is compatible with \(\Delta(c_1,d_1)\Delta(c_2,d_2)\). \(\Delta(1,1)\Delta(1,1) = \Delta(3,0)\); \(\Box(1,2)\Box(1,2) = \Box(4,3)\); \(\Box(1,2)\Box(2,1) = \Box(5,0)\).
Gallery¶
A standard (method 1) 8-frequency geodesic division of an icosahedron.
|
(4,3) GC operation on a cube, in canonical form. |
Standard (method 1) (5,3) subdivision of a tetrahedron. (Notice the
vertices in the upper left and lower right.) |
(5,3) subdivision of a tetrahedron, but using nslerp projection and
optimized to reduce variation in face area. |
Conway’s Game of Life on the surface of a (5,3) subdivided cube. Courtesy of Roger Kaufman. |
Repeatedly applying the (1,1) triangular subdivision using the flat
projection and not normalizing gives this interesting organic shape.
It resembles a viral capsid (moreso than geodesic polyhedra usually
resemble a viral capsid), or a strange fruit. Produced by
|
Improper geodesic polyhedron on the 3-dihedron. Note dangling faces.
|
GC operation on the 4-dihedron.
|
Quadrilateral nslerp2 preserves diagonal lines across the
face. |
An attractive floral pattern produced by
|
 A triangular balloon polyhedron. |
 A quadrilateral balloon polyhedron, which resembles the shape of a
peeled coconut. |
Indices and tables¶
References
| [Antiprism] | http://www.antiprism.com/ |
| [Polyhedronisme] | http://levskaya.github.io/polyhedronisme/ |
| [Brinkmann] | https://arxiv.org/abs/1705.02848 |
| [Deza2004] | M. Deza and M. Dutour. Goldberg-Coxeter constructions for 3- and 4-valent plane graphs. The Electronic Journal of Combinatorics, 11, 2004. #R20. |
| [Deza2015] | Michel-Marie Deza, Mathieu Dutour Sikirić, Mikhail Ivanovitch Shtogrin. (2015) Geometric Structure of Chemistry-Relevant Graphs: Zigzags and Central Circuits. Springer. pp 131-148. https://books.google.com/books?id=HLi4CQAAQBAJ&lpg=PA130&ots=ls1r5QkM51&dq=goldberg-coxeter&pg=PA130#v=onepage&q=goldberg-coxeter&f=false |
| [Altschuler] | Altschuler, E. L. et al. 1997. Possible Global Minimum Lattice Configurations for Thomson’s Problem of Charges on a Sphere. Phys. Rev. Lett. 78, 2681. [http://dx.doi.org/10.1103/PhysRevLett.78.2681 doi:10.1103/PhysRevLett.78.2681] http://www.mcs.anl.gov/~zippy/publications/thomson/thomsonPRL.html |
| [Goldberg] | Goldberg, M, 1937. A class of multi-symmetric polyhedra. Tohoku Mathematical Journal. |
| [Hart1997] | Hart, G, 1997. Calculating Canonical Polyhedra. Mathematica in Education and Research, Vol 6 No. 3, Summer 1997, pp. 5-10. |
| [Hart2012] | Hart, G, 2012. Goldberg Polyhedra. In Senechal, Marjorie. Shaping Space (2nd ed.). Springer. pp. 125–138. [http://dx.doi.org/10.1007/978-0-387-92714-5_9 doi:10.1007/978-0-387-92714-5_9]. |
| [Kenner] | Kenner, H, 1976. Geodesic Math and How to Use It. University of California Press. |
| [Schein] | Schein & Gayed, 2014. Fourth class of convex equilateral polyhedron with polyhedral symmetry related to fullerenes and viruses. PNAS, Early Edition [http://dx.doi.org/10.1073/pnas.1310939111 doi:10.1073/pnas.1310939111] |
| [Oosterom] | Van Oosterom, A & Strackee, J, 1983. The Solid Angle of a Plane Triangle. IEEE Trans. Biom. Eng. BME-30 (2): 125–126. [http://dx.doi.org/10.1109/TBME.1983.325207 doi:10.1109/TBME.1983.325207]. |
| [Folke] | Eriksson, Folke (1990). “On the measure of solid angles”. Math. Mag. 63 (3): 184–187. doi:10.2307/2691141. JSTOR 2691141. |
| [Caspar] | D.L.D. Caspar and A. Klug. Physical principles in the construction of regular viruses. In Cold Spring Harb Symp Quant Biol., volume 27, pages 1–24, 1962. |
| [Coxeter] | H.S.M. Coxeter. Virus macromolecules and geodesic domes. In J.C. Butcher, editor, A spectrum of mathematics, pages 98–107. Oxford University Press, 1971. |
| [Coxeter8] | H.S.M. Coxeter. Truncation. In Regular Polytopes, pages 145-164 3rd ed, Dover, 1973. |
| [Tarnai] | T. Tarnai, F. Kovacs, P.W. Fowler, and S.D. Guest. Wrapping the cube and other polyhedra. Proceedings of the Royal Society A, 468:2652–2666, 2012. |
| [StamLoop] | Jos Stam: Evaluation of Loop Subdivision Surfaces, Computer Graphics Proceedings ACM SIGGRAPH 1998, |
| [StamCatmull] | Stam, J. (1998). “Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values”. Proceedings of the 25th annual conference on Computer graphics and interactive techniques - SIGGRAPH ‘98. pp. 395–404. doi:10.1145/280814.280945. ISBN 0-89791-999-8. |
| [HartPropeller] | Hart, G. W. (2000) Sculpture Based on Propellorized Polyhedra. Proceedings of MOSAIC 2000, Seattle, Washington, August 2000, pp. 61-70. http://www.georgehart.com/propello/propello.html |
| [HartConway] | Hart, G. W. Conway Notation for Polyhedra. http://www.georgehart.com/virtual-polyhedra/conway_notation.html |