Calculating Canonical Polyhedra

A relaxation algorithm is presented to determine a canonical form for an arbitrary convex polyhedron.

by George W. Hart

One important role of high-level mathematical software such as Mathematica is that it easily allows for the testing of experimental algorithms. Here we explore a method of finding a canonical form of a polyhedron. The canonical form is a polyhedron topologically equivalent to an input polyhedron, but with all edges tangent to a unit sphere and with the center of gravity of the tangent points being the origin.  A variety of examples will illustrate this procedure.   One application is that the algorithm constructs a geometrically self-dual polyhedron given one which is only combinatorially self-dual.

A theorem of Schramm [Schramm 1992] states that for any given convex polyhedron or 3-connected planar graph, there is a topologically equivalent polyhedron with the following properties:

 1) each edge is tangent to the unit sphere,

 2) the center of gravity of the points of tangency is the origin.



The solution is unique, up to rotations and reflections of the sphere and so provides a kind of canonical representation.  Although we will not pursue it here, the theorem is actually more general in that it allows the sphere to be replaced with an arbitrary smooth convex body, e.g., an egg-shaped quasi-ellipsoid.  The theorem is closely related to the Koebe-Andreev-Thurston Circle Packing theorem for planar graphs.  See Ziegler [Ziegler 1995, p. 118] for discussion and other references. 



For example, if our input polyhedron is a geometrically distorted (but topologically unchanged) form of any Platonic or Archimedean solid, and we calculate its canonical form, our algorithm should output its undistorted form centered at the origin, with a midradius of 1, and with an arbitrary rotation.  E.g., given any parallelepiped as input, the canonical form output is a cube of edge (Sqrt[2]). For an arbitrary polyhedron, the canonical form is a way of illustrating its combinatorial or topological structure, which often lets one immediately see and understand its structure and symmetry.  



As an illustration, we will use the algorithm to find three geometrically self-dual polyhedra, given starting points which are only combinatorially self-dual.  It follows from the above that the dual polyhedron also shares properties (1) and (2).  Thus, for these three self-dual examples, we can make a compound of the polyhedron with its dual, i.e., with itself, in which both polyhedra have the same points of tangencies and at these points their edges cross each other's at right angles. 



The proofs cited above are of existence and not constructive, so we are interested in an algorithm for determining the canonical form of an input polyhedron. The algorithm proposed here operates by relaxation to iteratively move the vertices of the given polyhedron along a trajectory which converges at the canonical form. Although we begin and end with a polyhedron, during this relaxation, the object defined by the vertices is likely not to be a polyhedron geometrically.  Sets of points which belong to a face (combinatorially) are likely not to be coplanar. So we add a third condition for a solution, to the two above:



3) the faces are planar.

Algorithm

Our algorithm inputs a polyhedron and iteratively adjusts its vertices to slightly improve its conformance with the three conditions above. Within a couple of hundred iterations, it typically finds all the conditions are satisfied within a small tolerance, and stops. Three simple operations are all that is needed, corresponding to the three conditions:



    A) For each edge, the closest point to the origin is determined; call it p. If p lies at unit distance from the origin, condition (1) is satisfied for that edge. If not, a small fraction of p is added to the two vertices which define the edge, (in proportion to the sign and amount of the error) so that at the next iteration the edge will be closer to tangency with the unit sphere. 

    B) The center of gravity of all the points p is determined. If it is zero, condition (2) is satisfied.  If not, it is subtracted from all the vertices. 

    C) For each face, if the vertices lie in some plane in space, condition (3) is satisfied. If not, a plane which approximates it is computed. Each vertex of the face is then moved along a normal towards the plane. 



Iterating, it sometimes takes many steps for a correction at one part of the polyhedron to percolate its way around and equalize the conditions everywhere. On the examples below, between 50 and 100 iterations were sufficient to have all conditions satisfied within a tolerance of (10^(-8)). Although the individual steps are quite simple, this takes a minute or so in these examples and would take longer on more complex polyhedra. Refinement and optimization are left as future work.

Any canonicalization program must begin with a data structure for representing polyhedra positioned in space.  Our data structure is naturally of a two-part nature, consisting of the topological connectivity information (which remains constant) and a list of vertex coordinates (which is iteratively updated). We use v for the list of vertex coordinates, each a triple {x,y,z} of Cartesian components.  The topological information is given as a list f of faces, where each face is a list of vertex indices (in either clockwise or counter-clockwise order).  For example, a cube of edge 2 can be represented with the following v and f:

vCube = {{-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, {-1, 1, 1}, {1, -1, -1}, 
    {1, -1, 1}, {1, 1, -1}, {1, 1, 1}};

fCube = {{1, 2, 4, 3}, {5, 6, 8, 7}, {1, 2, 6, 5}, {3, 4, 8, 7}, 
    {1, 3, 7, 5}, {2, 4, 8, 6}};

A polyhedron in this v and f form can be displayed with the following view function, e.g. view[vCube, fCube].

view[v_, f_] := 
  Show[Graphics3D[Polygon /@ Map[v[[#1]] & , f, {2}]], PlotRange-> All,Boxed -> False]

As we will be manipulating edges, we write a function face2edges which creates a list of edges from a given face, and a function f2e which creates a list of all edges from the list of faces, f.  The use of Sort and Union eliminates duplicates, e.g., resulting in only twelve edges for the cube.

face2edges[face_] := MapThread[Sort[{#1, #2}] & , {face, RotateLeft[face]}];

f2e[f_] := Union @@ face2edges /@ f

eCube = f2e[fCube]

Given any line, the closest point on that line to the origin will be at distance 1 from the origin if the line is tangent to the unit sphere, i.e., an edge satisfying condition (1) above.  The following method finds the closest point on a line to the origin, when we know two points, P1 and P2, of the line.  Use a parameter alpha to represent an arbitrary point on the line as P1 + alpha L, where L=P2-P1 is a vector parallel to the line.  Observe that if p is the closest point, the vector p = P1 + alpha L is orthogonal to L.  Expressing this as a dot product being zero and solving for alpha gives alpha = - L.P1 / L.L, which is used in the following function.

closest[{P1_, P2_}] := 
      With[{L = P2 - P1}, P1 - ((L . P1 L)/L) . L]

The following statement finds a list of these points corresponding to the list of edges, which in this example gives us the cube's twelve edge midpoints.  The inner Map replaces indices into v with xyz coordinates.

closest /@ Map[vCube[[#1]] & , eCube, {2}]

For step (B) of the algorithm, we subtract from each vertex the mean of these points of tangency.

recenter[v_, e_] := 
  With[{centroid = Plus @@ closest /@ Map[v[[#1]] & , e, {2}]/Length[e]}, 
   (#1 - centroid & ) /@ v]

For step (A) we Scan through the edges, compute the distance of the tangent point from the origin, and make a correction to the two endpoints in proportion to how the distance differs from unity.  The parameter 0.5 is experimental; it is designed to move the points slowly, preventing large unstable dynamics.

tangentify[v_, e_] := 
  Module[{newV = v, t, c}, Scan[(t = closest[v[[#1]]]; 
        c = 0.5*(1 - Sqrt[t . t])*t; newV[[#1[[1]]]] += c; 
        newV[[#1[[2]]]] += c; ) & , e]; newV]

For step (C) we characterize an approximating face plane by its unit normal vector, n, and its distance from the origin. A fast approximation to the normal is just to sum the normals at each face vertex, obtained by taking the cross product of the two edge vectors meeting there.  The distance is obtained by projecting the centroid of the vertices on a unit vector in that direction.  (A least-squares fit of a plane to the points might perform better, but would be slower.)  Again, a stabilizing factor is employed, here 0.2.  Future work might optimize these factors.

unit[x_] := With[{mag2 = x . x}, If[mag2 != 0, x/Sqrt[mag2], x]]

cross[{ax_, ay_, az_}, {bx_, by_, bz_}] := 
  {ay*bz - az*by, az*bx - ax*bz, ax*by - ay*bx}

approxNormal[face_] := 
  unit[Plus @@ MapThread[unit[cross[#1 - #2, #2 - #3]] & , 
     {face, RotateLeft[face, 1], RotateLeft[face, 2]}]]

planarize[v_, f_] := 
  Module[{newV = v, faceXYZ, n, centroid}, 
   Scan[(faceXYZ = v[[#1]]; n = approxNormal[faceXYZ]; 
        centroid = Plus @@ faceXYZ/Length[faceXYZ]; 
        if[n . centroid < 0, n = -n]; 
        Scan[newV[[#1]] += 0.2*n . (centroid - v[[#1]])*n & , #1]; ) & , f]; 
    newV]

Assembling the above pieces into an iterative loop gives the following canonicalization function.  As a sanity test, we verify that the procedure operates on a cube of edge 2 to give a cube of edge (Sqrt[2]), i.e., with edges tangent to the unit sphere.

canonicalize[v_, f_] := 
  Module[{newV = N[v], e = f2e[f], oldV, maxChange}, 
   Do[oldV = newV; newV = tangentify[newV, e]; newV = recenter[newV, e]; 
      newV = planarize[newV, f]; maxChange = Max[Abs[oldV - newV]]; 
      If[maxChange < 10.^(-8), Break[]], {i, 80}]; 
    Print["Solved within ", maxChange]; newV]

canonicalize[vCube, fCube]

Examples

We begin with a polyhedron which John Conway has named a "hermaphrodite." It is composed of an n-gonal prism joined to an n-gonal pyramid, and so has an n-gon for a base, n rectangular sides attached to the base, and n triangular faces which meet at an apex.  Note that it is combinatorially self-dual.  Numbering the vertices from 1 to n around the base, then from n+1 to 2n around the juncture, and numbering the apex 2n+1, the following construction functions are straightforward.

vHermaphrodite[n_] := 
  N[Join[Table[{Cos[(2*Pi*i)/n], Sin[(2*Pi*i)/n], -1}, {i, n}], 
    Table[{Cos[(2*Pi*i)/n], Sin[(2*Pi*i)/n], 0}, {i, n}], {{0, 0, 1}}]]

fHermaphrodite[n_] := 
  Join[{Range[n], {1, n + 1, 2*n, n}, {n + 1, 2*n, 2*n + 1}}, 
   Table[{i, i + 1, n + i + 1, n + i}, {i, n - 1}], 
   Table[{n + i, n + i + 1, 2*n + 1}, {i, n - 1}]]

view[vHermaphrodite[7], fHermaphrodite[7]];

Figure 1. "Hermaphrodite" formed by joining a 7-sided prism and a 7-
  sided pyramid has 15 faces and 15 vertices.  The net of this polyhedron is 
    the same as the net of its dual, so it is combinatorially self-dual.

The above is not canonical.  We now canonicalize it and view the result.  As the vertices are adjusted to make the edges tangent to the unit sphere, the squares are changed to trapezoids.

fH7 = fHermaphrodite[7]; 
vH7c = canonicalize[vHermaphrodite[7], fH7]; 
view[vH7c, fH7];

Figure 2. "Hermaphrodite" now in canonical form is not only combinatorially self-
  dual, but also geometrically self-dual.

The canonical form is of special interest because it is geometrically self-dual.  Conveniently (because 7 is odd) we can position its dual by simple inversion through the origin, i.e., negating all the vertex coordinates.  With the original vertices numbered from 1 to 15 and the dual's vertices numbered from 16 to 30, a compound of the canonical form with its dual is produced as follows.  Note that every edge crosses an edge of the dual at right angles.  The edge intersections are the points where the edges are tangent to the unit sphere.

view[Join[vH7c, -vH7c], Join[fH7, 15 + fH7]];

Figure 3. Combining the canonical hermaphrodite with itself properly rotated gives a 
  compound.  The duality of the two is evidenced by the fact that the edges 
  cross at right angles.  The crossing points all lie on the unit sphere.

A second example Conway has called the "antihermaphrodite," composed of half an n-gonal antiprism joined to half an n-gonal trapezohedron (the dual of an antiprism).  It has an n-gon for a base, n triangular sides attached to the base, and n kite-shaped faces which meet at an apex.  It is also self-dual, so in its canonical form, it makes an attractive compound with itself.  Numbering the vertices as above, the following code is analogous. We do not view this before canocalization however, as this funcion only provides a combinatorial polyhedron, in which the 4-sided faces are not planar.  After canonicalization we have a true geometric polyhedron with planar faces.

vAntiHermaphrodite[n_] := 
  N[Join[Table[{Cos[2*Pi*i/n], Sin[2*Pi*i/n], -1}, {i, n}], 
    Table[{Cos[2*Pi*(i + 0.5)/n], Sin[2*Pi*(i + 0.5)/n], 0}, {i, n}], 
    {{0, 0, 1}}]]

fAntiHermaphrodite[n_] := 
  Join[{Range[n], {1, 2*n, n}, {1, n + 1, 2*n + 1, 2*n}}, 
   Table[{i, i + 1, n + i}, {i, n - 1}], 
   Table[{i, n + i, 2*n + 1, n + i - 1}, {i, 2, n}]]

fA7 = fAntiHermaphrodite[7]; vA7c = canonicalize[vAntiHermaphrodite[7], fA7]; 
  view[vA7c, fA7];

Figure 4. "Antihermaphrodite" in canonical form also has 15 sides and 15 
  vertices.  It can be seen as half an antiprism combined with half a 
  trapezohedron.

Again this is geometrically self-dual, and we may compound the result with its dual, and see that every edge crosses an edge of the dual at right angles.  In this case, we position the dual by negating the z-component of each vertex.

view[Join[vA7c, ({#1[[1]], #1[[2]], -#1[[3]]} & ) /@ vA7c], 
  Join[fA7, 15 + fA7]]

Figure 5.  Combining the antihermaphrodite in canonical form with itself properly 
  rotated gives a compound.  Again the edges cross at right angles and the 
  crossing points all lie on the unit sphere, illustrating self-duality

As a final example, we create a "tetrahedrally stellated icosahedron" (TSI) by choosing four faces of an icosahedron in the face planes of a tetrahedron and covering them with low triangular pyramids.  If the height of the pyramid is chosen properly,  the pyramid faces become an extensions of the adjacent icosahedron faces.  The resulting polyhedron has twelve kite-shaped faces and four triangular faces, and is an interesting model of chiral tetrahedral symmetry, i.e., it has four 3-fold axes and three 2-fold axes, but no mirror planes.   Examining its net shows that the TSI is combinatorially self-dual.  We construct it from the vertices of an icosahedron and a surrounding tetrahedron, manually writing the vertex indices from a pencil sketch:

vTSI = With[{t = 1.618, d = 1.1708}, 
    {{-1, t, 0}, {1, t, 0}, {-1, -t, 0}, {1, -t, 0}, {0, 1, -t}, {0, -1, -t}, 
     {0, 1, t}, {0, -1, t}, {-t, 0, -1}, {-t, 0, 1}, {t, 0, -1}, {t, 0, 1}, 
     {d, d, d}, {d, -d, -d}, {-d, d, -d}, {-d, -d, d}}];

fTSI = {{1, 2, 5, 15}, {1, 7, 13, 2}, {3, 4, 14, 6}, {3, 4, 8, 16}, 
    {5, 6, 9, 15}, {5, 6, 14, 11}, {7, 8, 16, 10}, {7, 8, 12, 13}, 
    {9, 10, 1, 15}, {9, 10, 16, 3}, {11, 12, 13, 2}, {11, 12, 4, 14}, 
    {1, 7, 10}, {2, 5, 11}, {3, 6, 9}, {4, 8, 12}};

view[vTSI, fTSI];

Figure 6. "Tetrahedrally stellated icosahedron" (TSI) is based on an icosahedron with
   four faces covered with low triangular pyramids that blend in with the 
  adjoining icosahedral faces.  It has 16 faces and 16 vertices.

Next, we find its canonical form, which has a different shape for the kite faces, and is geometrically self-dual.

vTSIc = canonicalize[vTSI, fTSI]; 
view[vTSIc, fTSI];

Figure 7. Canonicalized TSI output by the algorithm.  The twelve quadrilateral faces 
  have a different shape than in Figure 6.

To position this in a compound with its dual, we have to rotate it 90 degrees about one of its 2-fold axes.  By the above construction, these are the x, y, and z axes, and we choose the y-axis here:

view[Join[vTSIc, ({-#1[[3]], #1[[2]], #1[[1]]} & ) /@ vTSIc], 
  Join[fTSI, 16 + fTSI]]

Figure 8.  Combining the canonical TSI with a rotated version of itself gives a 
  compound, again illustrating self-duality.

Conclusions and future work

A polyhedron canonicalization algorithm has been presented and demonstrated with a range of examples.   See [Hart, 1997] for additional examples.  Empirically, the algorithm has always  been observed to converge, although we have no proof that it should do so generally.  One application demonstrated is to construct geometrically self-dual polyhedra which form attractive compounds with themselves.  Three-dimensional virtual reality models of these and other examples are available on the electronic version of this journal.



There are many directions for extending this work, e.g., faster algorithms, provably convergent algorithms, or tangency to an arbitrary, nonspherical, smooth, convex solid. It seems likely that a "primal-dual" algorithm which simultaneously relaxes both the original polyhedron and its dual should be feasible and have certain advantages.

Acknowledgment

John Conway told me of this theorem, the need for a construction algorithm, and the self-dual hermaphrodite and antihermaphrodite polyhedra.

References

Coxeter, H.S.M. , Regular Polytopes, Macmillan, 1963, (Dover reprint, 1973).

Ziegler, G. M., Lectures on Polytopes, Springer-Verlag, 1995.

Schramm, O., "How to Cage an Egg," Inventiones Math. 107, 1992, pp. 543-560.

Hart, G. W., Virtual Polyhedra, 1997, http://www.li.net/~george/virtual-polyhedra/vp.html

ABOUT THE AUTHOR

George W. Hart
george@li.net
http://www.li.net/~george
Department of Computer Science
Hofstra University
Hempstead, NY 11550