Truncated icosahedron

truncated_icosahedron.m

clear all;

# The dot products between the symmetry plane normals.
plane_dots = [
          1   -cos(pi/3)           0
  -cos(pi/3)           1  -cos(pi/5)
          0   -cos(pi/5)           1
];

# The dot products of the seed vertices with the symmetry plane normals.
seed_dots = [ 1 1 0; 1 0 0; 0 0 1 ].';

# [ n1x n1y n1z ] [ n1x n2x n3x ]
# [ n2x n2y n2z ] [ n1y n2y n3y ] = plane_dots
# [ n3x n3y n3z ] [ n1z n2z n3z ]

# [ n1x   0   0 ] [ n1x n2x n3x ]
# [ n2x n2y   0 ] [   0 n2y n3y ] = plane_dots
# [ n3x n3y n3z ] [   0   0 n3z ]

# n1x^2 = 1
# n1x n2x = -cos(pi/3)
# n1x n3x = 0
# n2x^2 + n2y^2 = 1
# n2x n3x + n2y n3y = -cos(pi/5)
# n3x^2 + n3y^2 + n3z^2 = 1

# n1x = 1
# n2x = -cos(pi/3)
# n2y = sin(pi/3)
# n3x = 0
# n3y = -cos(pi/5) / sin(pi/3)
# n3z = sqrt(1 - (-cos(pi/5) / sin(pi/3))^2)

# plane_normals * plane_normals.' = plane_dots
# R' * R = A
# R = chol(A)
plane_normals = chol(plane_dots).';
plane_normals

# [ n1x   0   0 ] [ sx ]   [ x ]
# [ n2x n2y   0 ] [ sy ] = [ y ]
# [ n3x n3y n3z ] [ sz ]   [ z ]

# [          1                     0                                 0 | x ]
# [ -cos(pi/3)             sin(pi/3)                                 0 | y ]
# [          0  -cos(pi/5)/sin(pi/3)  sqrt(1-(-cos(pi/5)/sin(pi/3))^2) | z ]

# [ 1  0  0 | x                                                                        ]
# [ 0  1  0 | (cos(pi/3)x+y)/sin(pi/3)                                                 ]
# [ 0  0  1 | (cos(pi/5)(cos(pi/3)x+y)/sin(pi/3)^2+z)/sqrt(1-(-cos(pi/5)/sin(pi/3))^2) ]

# plane_normals * seed_vertex = seed_dots
seed_vertices = inv(plane_normals) * seed_dots;
seed_vertices = seed_vertices ./ norm(seed_vertices, "columns");
seed_vertices

vertices = [ seed_vertices ];
symmetry_queue = vertices;

while (!isempty(symmetry_queue))
  vertex = symmetry_queue(:, end);
  symmetry_queue(:, end) = [];

  projected_to_normals = plane_normals .* (plane_normals * vertex);
  mirrored_vertices = (vertex.' .- 2 * projected_to_normals).';

  for new_vertex = mirrored_vertices
    vertex_deltas = vertices .- new_vertex;
    vertex_delta_sq_mags = dot(vertex_deltas, vertex_deltas, 1);
    if min(vertex_delta_sq_mags) >= 1e-6
      # The vertex is not in the output yet.
      vertices(:, end+1) = new_vertex;
      symmetry_queue(:, end+1) = new_vertex;
    endif
  endfor
endwhile

newplot();
hold on;
#scatter3(vertices(1,:), vertices(2,:), vertices(3,:));
scatter3(seed_vertices(1,:), seed_vertices(2,:), seed_vertices(3,:));
k = convhulln(vertices.');
trisurf(k, vertices(1,:), vertices(2,:), vertices(3,:),
  "FaceColor", "cyan", "EdgeColor", "none", "FaceLighting", "Gouraud");
hold off;

axis([-1 1 -1 1 -1 1], "square");
rotate3d on;
shading interp;
lighting flat;
light();