Here is a slight shortening of Clayton's code:

DynamicModule[{p, b, verts = Normalize /@ N[PolyhedronData["Cube", "VertexCoordinates"]],
               t = ArcCos[1/3]/2, cols = RGBColor /@ {"#EF6C35", "#2BB3C0", "#161C2E"}}, 
              Manipulate[p = verts.RotationMatrix[-?, {1, 0, 0}]; b = NullSpace[{#}] & /@ p;
                         Graphics[{EdgeForm[], Table[{cols[[Quotient[i, 4, -3]]], 
                                   Polygon[Table[Most[#]/(1 - Last[#]) &[
                                                 Cos[t] p[[i]] + Sin[t] ({Cos[s], Sin[s]}.b[[i]])],
                                                 {s, 0, 2 ?, 2 ?/200}]]}, {i, Length[verts]}]}, 
                                  PlotRange -> Sqrt[6], Background -> cols[[-1]], ImageSize -> 540],
                         {?, 0, ?/2}]]

Notes:

1. Dot products are very efficient in Mathematica, thus the new direct expression for p (accounting for the need to transpose) is much faster than mapping a dot product all over a list.
2. Using inexact coordinates in this case, apart from the easier internal handling, allows NullSpace[] to use the singular value decomposition (SVD) internally, thus obviating the need for an extra Orthogonalize[].