Griddip (truncated icosidodecahedron) with rgl (R)

griddip.R

library(cxhull)
library(rgl)

rotation4Dplane <- function(axis1, axis2, theta, vector){
  # axes must be normalized
  vx <- sum(vector*axis1)
  vy <- sum(vector*axis2)
  cos_theta <- cos(theta)
  sin_theta <- sin(theta)
  coef1 <- vx*cos_theta - vy*sin_theta
  coef2 <- vy*cos_theta + vx*sin_theta
  pvector <- vx*axis1 + vy*axis2 
  coef1*axis1 + coef2*axis2 + (vector-pvector)
}


#~~~~ truncated icosidodecahedron with edge length 2*phi-2 ~~~~#
phi <- (1 + sqrt(5)) / 2
vs0 <- rbind(
  c(   1/phi, 1/phi,    3+phi ),
  c(   2/phi,   phi,  1+2*phi ),
  c(   1/phi, phi^2, -1+3*phi ),
  c( 2*phi-1,     2,    2+phi ),
  c(     phi,     3,    2*phi )
)
# changes of sign
vs1 <- gyro::changesOfSign(vs0)
# even permutations
vs2 <- rbind(
  vs1, vs1[, c(3L, 1L, 2L)], vs1[, c(2L, 3L, 1L)]
)

# Cartesian product with a segment
vertices <- rbind(
  cbind(vs2, -phi+1), 
  cbind(vs2, phi-1)
)
summary(apply(vertices, 1L, crossprod)) # all equal
r <- sqrt(c(crossprod(vertices[1L, ])))

# hull
library(cxhull)
hull <- cxhull(vertices, triangulate = FALSE)
cells <- hull[["facets"]]
ridges <- hull[["ridges"]]

# cubic cells
cubicCells <- Filter(function(f) length(f[["vertices"]]) == 8, cells)

# squares and their edges
goodRidges <- do.call(c, lapply(cubicCells, `[[`, "ridges"))
polygonize <- function(edges){ # function which orders the vertices
  nedges <- nrow(edges)
  indices <- edges[1L, ]
  i <- indices[2L]
  edges <- edges[-1L, ]
  for(. in 1L:(nedges-2L)){
    j <- which(apply(edges, 1L, function(e) i %in% e))
    i <- edges[j, ][which(edges[j, ] != i)]
    indices <- c(indices, i)
    edges <- edges[-j, ]
  }
  indices
}
squares <- t(vapply(
  goodRidges, 
  function(r) polygonize(ridges[[r]][["edges"]]), integer(4L)
))
edges <- do.call(
  rbind,
  apply(squares, 1L, function(x){
    cbind(x, c(x[-1L], x[1L]))
  }, simplify = FALSE)
)
edges <- t(apply(edges, 1L, sort))
edges <- edges[!duplicated(edges), ]

# Stereographic-like projection
sproj <- function(v, r){
  acos(v[4]/r) / (r^4 - v[4L]^4)^(1/4) * v[1L:3L]
}

# rotated and projected vertices
verts3D <- function(xi){
  rpoints <- apply(vertices, 1L, function(v){
    rotation4Dplane(c(1, 0, 0, 0), c(0, 1, 0, 0), xi, v)
  })
  t(apply(rpoints, 2L, function(point) sproj(point, r)))
}


# rgl ####
xi_ <- seq(0, 2*pi, length.out = 121L)[-1L]
open3d(windowRect = c(50, 50, 562, 562), zoom = 0.65)
bg3d(rgb(54, 57, 64, maxColorValue = 255))
for(j in seq_along(xi_)){
  spheres3d(0, 0, 0, radius = 2, alpha = 0) # bounding sphere
  points <- verts3D(xi_[j])
  for(i in 1L:nrow(squares)){
    quads3d(points[squares[i, ], ], color = "darkorange")
  }
  for(i in 1L:nrow(edges)){
    edge <- edges[i, ]
    shade3d(cylinder3d(
      rbind(
        points[edge[1L], ],
        points[edge[2L], ]
      ), 
      radius = 0.02, sides = 60L
    ), color = "gold")
  }
  spheres3d(points, radius = 0.03, color = "gold")
  rgl.snapshot(sprintf("zzz_griddip%03d.png", j))
  clear3d()
}

# make GIF ####
library(gifski)
pngs <- list.files(pattern = "^zzz_griddip")
gifski(
  pngs,
  gif_file = "griddip.gif",
  width = 512, height = 512,
  delay = 1/10,
  loop = TRUE
)

file.remove(pngs)

griddip.gif