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.