Skip to content

Instantly share code, notes, and snippets.

@stla

stla/runcinatedTesseract.R

Last active August 8, 2023 10:21
Show Gist options
  • Save stla/61206316a3bc5a38c7268cec7ff3ca20 to your computer and use it in GitHub Desktop.
Save stla/61206316a3bc5a38c7268cec7ff3ca20 to your computer and use it in GitHub Desktop.
Download ZIP
Runcinated tesseract (tetrahedra only) with R
Raw
runcinatedTesseract.R
library(rgl)
library(Rvcg)
# vertices ####
nvertices <- 64L
vertices <-
rbind(
c(-1, -1, -1, -1 - sqrt(2)),
c(1, -1, -1, -1 - sqrt(2)),
c(-1, 1, -1, -1 - sqrt(2)),
c(1, 1, -1, -1 - sqrt(2)),
c(-1, -1, 1, -1 - sqrt(2)),
c(1, -1, 1, -1 - sqrt(2)),
c(-1, 1, 1, -1 - sqrt(2)),
c(1, 1, 1, -1 - sqrt(2)),
c(-1, -1, -1, 1 + sqrt(2)),
c(1, -1, -1, 1 + sqrt(2)),
c(-1, 1, -1, 1 + sqrt(2)),
c(1, 1, -1, 1 + sqrt(2)),
c(-1, -1, 1, 1 + sqrt(2)),
c(1, -1, 1, 1 + sqrt(2)),
c(-1, 1, 1, 1 + sqrt(2)),
c(1, 1, 1, 1 + sqrt(2)),
c(-1, -1, -1 - sqrt(2), -1),
c(1, -1, -1 - sqrt(2), -1),
c(-1, 1, -1 - sqrt(2), -1),
c(1, 1, -1 - sqrt(2), -1),
c(-1, -1, 1 + sqrt(2), -1),
c(1, -1, 1 + sqrt(2), -1),
c(-1, 1, 1 + sqrt(2), -1),
c(1, 1, 1 + sqrt(2), -1),
c(-1, -1, -1 - sqrt(2), 1),
c(1, -1, -1 - sqrt(2), 1),
c(-1, 1, -1 - sqrt(2), 1),
c(1, 1, -1 - sqrt(2), 1),
c(-1, -1, 1 + sqrt(2), 1),
c(1, -1, 1 + sqrt(2), 1),
c(-1, 1, 1 + sqrt(2), 1),
c(1, 1, 1 + sqrt(2), 1),
c(-1, -1 - sqrt(2), -1, -1),
c(1, -1 - sqrt(2), -1, -1),
c(-1, 1 + sqrt(2), -1, -1),
c(1, 1 + sqrt(2), -1, -1),
c(-1, -1 - sqrt(2), 1, -1),
c(1, -1 - sqrt(2), 1, -1),
c(-1, 1 + sqrt(2), 1, -1),
c(1, 1 + sqrt(2), 1, -1),
c(-1, -1 - sqrt(2), -1, 1),
c(1, -1 - sqrt(2), -1, 1),
c(-1, 1 + sqrt(2), -1, 1),
c(1, 1 + sqrt(2), -1, 1),
c(-1, -1 - sqrt(2), 1, 1),
c(1, -1 - sqrt(2), 1, 1),
c(-1, 1 + sqrt(2), 1, 1),
c(1, 1 + sqrt(2), 1, 1),
c(-1 - sqrt(2), -1, -1, -1),
c(1 + sqrt(2), -1, -1, -1),
c(-1 - sqrt(2), 1, -1, -1),
c(1 + sqrt(2), 1, -1, -1),
c(-1 - sqrt(2), -1, 1, -1),
c(1 + sqrt(2), -1, 1, -1),
c(-1 - sqrt(2), 1, 1, -1),
c(1 + sqrt(2), 1, 1, -1),
c(-1 - sqrt(2), -1, -1, 1),
c(1 + sqrt(2), -1, -1, 1),
c(-1 - sqrt(2), 1, -1, 1),
c(1 + sqrt(2), 1, -1, 1),
c(-1 - sqrt(2), -1, 1, 1),
c(1 + sqrt(2), -1, 1, 1),
c(-1 - sqrt(2), 1, 1, 1),
c(1 + sqrt(2), 1, 1, 1)
)
# tetrahedra ####
tetrahedra <-
rbind(
c(38, 22, 6, 54),
c(30, 14, 46, 62),
c(63, 31, 47, 15),
c(39, 55, 7, 23),
c(26, 10, 42, 58),
c(32, 16, 0, 48),
c(19, 35, 51, 3),
c(44, 12, 28, 60),
c(27, 59, 43, 11),
c(18, 2, 50, 34),
c(45, 29, 61, 13),
c(33, 17, 1, 49),
c(41, 25, 9, 57),
c(37, 5, 53, 21),
c(40, 24, 56, 8),
c(36, 4, 52, 20)
) + 1
# tetrahedra faces (triangles) ####
faces <- matrix(NA_real_, nrow = 64L, ncol = 3L)
for(i in 1:16) {
faces[4*(i-1) + 1, ] <- tetrahedra[i, c(0, 1, 2)+1]
faces[4*(i-1) + 2, ] <- tetrahedra[i, c(0, 1, 3)+1]
faces[4*(i-1) + 3, ] <- tetrahedra[i, c(0, 2, 3)+1]
faces[4*(i-1) + 4, ] <- tetrahedra[i, c(1, 2, 3)+1]
}
# edges ####
edges <- matrix(NA_real_, nrow = 96L, ncol = 2L)
for(i in 1:16) {
edges[6*(i-1) + 1, ] <- tetrahedra[i, c(0, 1)+1]
edges[6*(i-1) + 2, ] <- tetrahedra[i, c(0, 2)+1]
edges[6*(i-1) + 3, ] <- tetrahedra[i, c(0, 3)+1]
edges[6*(i-1) + 4, ] <- tetrahedra[i, c(1, 2)+1]
edges[6*(i-1) + 5, ] <- tetrahedra[i, c(1, 3)+1]
edges[6*(i-1) + 6, ] <- tetrahedra[i, c(2, 3)+1]
}
# 4D rotation ####
rotate4d <- function(alpha, beta, xi, vec){
a <- cos(xi)
b <- sin(alpha) * cos(beta) * sin(xi)
c <- sin(alpha) * sin(beta) * sin(xi)
d <- cos(alpha) * sin(xi)
x <- vec[1L]; y <- vec[2L]; z <- vec[3L]; w <- vec[4L]
c(
a*x - b*y - c*z - d*w,
a*y + b*x + c*w - d*z,
a*z - b*w + c*x + d*y,
a*w + b*z - c*y + d*x
)
}
# stereographic projection ####
stereographicProjection <- function(q, r2) {
acos(q[4]/sqrt(r2))/pi/sqrt(r2-q[4]*q[4]) * q[1:3]
}
# rotated and projected vertices ####
Vertices3 <- function(theta, phi, xi, r2) {
out <- matrix(NA_real_, nrow = nvertices, ncol = 3L)
for(i in 1:nvertices) {
out[i, ] <- stereographicProjection(
rotate4d(theta, phi, xi, vertices[i, ]),
r2
)
}
out
}
# spherical segment ####
sphericalSegment <- function(P, Q, r, n) {
out <- matrix(NA_real_, nrow = n+1, ncol = 4)
for(i in 0:n) {
pt <- P + (i/n)*(Q-P)
out[i+1, ] <- r / sqrt(c(crossprod(pt))) * pt
}
out
}
# draw an edge ####
Edge <- function(verts, v1, v2, theta, phi, xi, r2) {
P <- verts[v1, ]
Q <- verts[v2, ]
PQ <- sphericalSegment(P, Q, sqrt(r2), 100)
ctr <- matrix(NA_real_, nrow = 101L, ncol = 3L)
radius <- numeric(101L)
for(k in 1:101) {
O <- ctr[k, ] <- stereographicProjection(
rotate4d(theta, phi, xi, PQ[k, ]),
r2
)
radius[k] <- log1p(sqrt(c(crossprod(O)))/4)/5
}
cyl <- cylinder3d(ctr, radius, sides = 30)
shade3d(cyl, color = "yellow")
}
# stereographic subdivision ####
midpoint4 <- function(A, B, r) {
xmid <- (A[1] + B[1]) / 2
ymid <- (A[2] + B[2]) / 2
zmid <- (A[3] + B[3]) / 2
tmid <- (A[4] + B[4]) / 2
mids <- c(xmid, ymid, zmid, tmid)
lg <- sqrt(c(crossprod(mids))) / r
mids / lg
}
subdiv0 <- function(A, B, C, r) {
mAB <- midpoint4(A, B, r)
mAC <- midpoint4(A, C, r)
mBC <- midpoint4(B, C, r)
trgl1 <- rbind(A, mAB, mAC)
trgl2 <- rbind(B, mBC, mAB)
trgl3 <- rbind(C, mAC, mBC)
trgl4 <- rbind(mAB, mAC, mBC)
list(trgl1, trgl2, trgl3, trgl4)
}
subdiv <- function(A, B, C, r, depth) {
if(depth == 1) {
out <- subdiv0(A, B, C, r)
} else {
triangles <- subdiv(A, B, C, r, depth-1)
out <- vector("list", 4^(depth-1))
for(i in 1:(4^(depth-1))) {
trgl <- triangles[[i]]
trgls <- subdiv0(trgl[1, ], trgl[2, ], trgl[3,], r)
out[[4*(i-1)+1]] = trgls[[1]]
out[[4*(i-1)+2]] = trgls[[2]]
out[[4*(i-1)+3]] = trgls[[3]]
out[[4*(i-1)+4]] = trgls[[4]]
}
}
out
}
#############################
#### draw the polychoron ####
#############################
theta <- pi/2
phi <- 0
xi <- 0
r2 <- 3 + (1+sqrt(2))^2
r <- sqrt(r2)
depth <- 4
vs <- Vertices3(theta, phi, xi, r2)
stereoTriangles <- vector("list", 64)
for(i in 1:64) {
triangles4 <- subdiv(
vertices[faces[i, 1], ],
vertices[faces[i, 2], ],
vertices[faces[i, 3], ],
r, depth
)
stereoTriangles[[i]] <- vector("list", length(triangles4))
for(j in 1:length(triangles4)) {
trgl4 <- triangles4[[j]]
stereoTriangles[[i]][[j]] <-
rbind(
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[1, ]), r2
),
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[2, ]), r2
),
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[3, ]), r2
)
)
}
}
# plot ####
open3d(windowRect = 50 + c(0, 0, 512, 512))
view3d(0, 0, zoom = 0.9)
# draw triangles
for(i in 1:64) {
trgl <- stereoTriangles[[i]]
vrtcs <- t(do.call(rbind, trgl))
indices <- matrix(1L:(3L*length(trgl)), nrow = 3L)
mesh <- tmesh3d(
vertices = vrtcs,
indices = indices
)
mesh <- Rvcg::vcgClean(mesh, sel = c(0, 7), silent = TRUE)
shade3d(mesh, color = "red")
}
# draw edges
for(i in 1:96) {
Edge(vertices, edges[i, 1], edges[i, 2], theta, phi, xi, r2)
}
# draw vertices
for(i in 1:64) {
spheres3d(
rbind(vs[i,]),
radius = log1p(sqrt(c(crossprod(vs[i, ])))/4)/4,
color = "yellow"
)
}
# animation ####
M <- par3d("userMatrix")
movie3d(
par3dinterp(
time = seq(0, 1, len = 9),
userMatrix = list(
M,
rotate3d(M, pi, 1, 0, 0),
rotate3d(M, pi, 1, 1, 0),
rotate3d(M, pi, 1, 1, 1),
rotate3d(M, pi, 0, 1, 1),
rotate3d(M, pi, 0, 1, 0),
rotate3d(M, pi, 1, 0, 1),
rotate3d(M, pi, 0, 0, 1),
M
)
),
fps = 120,
duration = 1,
dir = ".",
frames = "zzzpic",
convert = FALSE,
clean = FALSE,
webshot = FALSE
)
pngs <- list.files(".", pattern = "^zzzpic", full.names = TRUE)
library(gifski)
gifski(pngs, "runcinatedTesseract_tetrahedra.gif",
width = 512, height = 512, delay = 1/9)
file.remove(pngs)
Raw
runcinatedTesseract2.R
library(rgl)
library(Rvcg)
# vertices ####
nvertices <- 64L
vertices <-
rbind(
c(-1, -1, -1, -1 - sqrt(2)),
c(1, -1, -1, -1 - sqrt(2)),
c(-1, 1, -1, -1 - sqrt(2)),
c(1, 1, -1, -1 - sqrt(2)),
c(-1, -1, 1, -1 - sqrt(2)),
c(1, -1, 1, -1 - sqrt(2)),
c(-1, 1, 1, -1 - sqrt(2)),
c(1, 1, 1, -1 - sqrt(2)),
c(-1, -1, -1, 1 + sqrt(2)),
c(1, -1, -1, 1 + sqrt(2)),
c(-1, 1, -1, 1 + sqrt(2)),
c(1, 1, -1, 1 + sqrt(2)),
c(-1, -1, 1, 1 + sqrt(2)),
c(1, -1, 1, 1 + sqrt(2)),
c(-1, 1, 1, 1 + sqrt(2)),
c(1, 1, 1, 1 + sqrt(2)),
c(-1, -1, -1 - sqrt(2), -1),
c(1, -1, -1 - sqrt(2), -1),
c(-1, 1, -1 - sqrt(2), -1),
c(1, 1, -1 - sqrt(2), -1),
c(-1, -1, 1 + sqrt(2), -1),
c(1, -1, 1 + sqrt(2), -1),
c(-1, 1, 1 + sqrt(2), -1),
c(1, 1, 1 + sqrt(2), -1),
c(-1, -1, -1 - sqrt(2), 1),
c(1, -1, -1 - sqrt(2), 1),
c(-1, 1, -1 - sqrt(2), 1),
c(1, 1, -1 - sqrt(2), 1),
c(-1, -1, 1 + sqrt(2), 1),
c(1, -1, 1 + sqrt(2), 1),
c(-1, 1, 1 + sqrt(2), 1),
c(1, 1, 1 + sqrt(2), 1),
c(-1, -1 - sqrt(2), -1, -1),
c(1, -1 - sqrt(2), -1, -1),
c(-1, 1 + sqrt(2), -1, -1),
c(1, 1 + sqrt(2), -1, -1),
c(-1, -1 - sqrt(2), 1, -1),
c(1, -1 - sqrt(2), 1, -1),
c(-1, 1 + sqrt(2), 1, -1),
c(1, 1 + sqrt(2), 1, -1),
c(-1, -1 - sqrt(2), -1, 1),
c(1, -1 - sqrt(2), -1, 1),
c(-1, 1 + sqrt(2), -1, 1),
c(1, 1 + sqrt(2), -1, 1),
c(-1, -1 - sqrt(2), 1, 1),
c(1, -1 - sqrt(2), 1, 1),
c(-1, 1 + sqrt(2), 1, 1),
c(1, 1 + sqrt(2), 1, 1),
c(-1 - sqrt(2), -1, -1, -1),
c(1 + sqrt(2), -1, -1, -1),
c(-1 - sqrt(2), 1, -1, -1),
c(1 + sqrt(2), 1, -1, -1),
c(-1 - sqrt(2), -1, 1, -1),
c(1 + sqrt(2), -1, 1, -1),
c(-1 - sqrt(2), 1, 1, -1),
c(1 + sqrt(2), 1, 1, -1),
c(-1 - sqrt(2), -1, -1, 1),
c(1 + sqrt(2), -1, -1, 1),
c(-1 - sqrt(2), 1, -1, 1),
c(1 + sqrt(2), 1, -1, 1),
c(-1 - sqrt(2), -1, 1, 1),
c(1 + sqrt(2), -1, 1, 1),
c(-1 - sqrt(2), 1, 1, 1),
c(1 + sqrt(2), 1, 1, 1)
)
# tetrahedra ####
tetrahedra <-
rbind(
c(38, 22, 6, 54),
c(30, 14, 46, 62),
c(63, 31, 47, 15),
c(39, 55, 7, 23),
c(26, 10, 42, 58),
c(32, 16, 0, 48),
c(19, 35, 51, 3),
c(44, 12, 28, 60),
c(27, 59, 43, 11),
c(18, 2, 50, 34),
c(45, 29, 61, 13),
c(33, 17, 1, 49),
c(41, 25, 9, 57),
c(37, 5, 53, 21),
c(40, 24, 56, 8),
c(36, 4, 52, 20)
) + 1
# tetrahedra faces (triangles) ####
faces <- matrix(NA_real_, nrow = 64L, ncol = 3L)
for(i in 1:16) {
faces[4*(i-1) + 1, ] <- tetrahedra[i, c(0, 1, 2)+1]
faces[4*(i-1) + 2, ] <- tetrahedra[i, c(0, 1, 3)+1]
faces[4*(i-1) + 3, ] <- tetrahedra[i, c(0, 2, 3)+1]
faces[4*(i-1) + 4, ] <- tetrahedra[i, c(1, 2, 3)+1]
}
# edges ####
edges <- matrix(NA_real_, nrow = 96L, ncol = 2L)
for(i in 1:16) {
edges[6*(i-1) + 1, ] <- tetrahedra[i, c(0, 1)+1]
edges[6*(i-1) + 2, ] <- tetrahedra[i, c(0, 2)+1]
edges[6*(i-1) + 3, ] <- tetrahedra[i, c(0, 3)+1]
edges[6*(i-1) + 4, ] <- tetrahedra[i, c(1, 2)+1]
edges[6*(i-1) + 5, ] <- tetrahedra[i, c(1, 3)+1]
edges[6*(i-1) + 6, ] <- tetrahedra[i, c(2, 3)+1]
}
# 4D rotation ####
rotate4d <- function(alpha, beta, xi, vec){
a <- cos(xi)
b <- sin(alpha) * cos(beta) * sin(xi)
c <- sin(alpha) * sin(beta) * sin(xi)
d <- cos(alpha) * sin(xi)
x <- vec[1L]; y <- vec[2L]; z <- vec[3L]; w <- vec[4L]
c(
a*x - b*y - c*z - d*w,
a*y + b*x + c*w - d*z,
a*z - b*w + c*x + d*y,
a*w + b*z - c*y + d*x
)
}
# stereographic projection ####
stereographicProjection <- function(q, r2) {
acos(q[4]/sqrt(r2))/pi/sqrt(r2-q[4]*q[4]) * q[1:3]
}
# rotated and projected vertices ####
Vertices3 <- function(theta, phi, xi, r2) {
out <- matrix(NA_real_, nrow = nvertices, ncol = 3L)
for(i in 1:nvertices) {
out[i, ] <- stereographicProjection(
rotate4d(theta, phi, xi, vertices[i, ]),
r2
)
}
out
}
# spherical segment ####
sphericalSegment <- function(P, Q, r, n) {
out <- matrix(NA_real_, nrow = n+1, ncol = 4)
for(i in 0:n) {
pt <- P + (i/n)*(Q-P)
out[i+1, ] <- r / sqrt(c(crossprod(pt))) * pt
}
out
}
# draw an edge ####
Edge <- function(verts, v1, v2, theta, phi, xi, r2) {
P <- verts[v1, ]
Q <- verts[v2, ]
PQ <- sphericalSegment(P, Q, sqrt(r2), 100)
ctr <- matrix(NA_real_, nrow = 101L, ncol = 3L)
radius <- numeric(101L)
for(k in 1:101) {
O <- ctr[k, ] <- stereographicProjection(
rotate4d(theta, phi, xi, PQ[k, ]),
r2
)
radius[k] <- log1p(sqrt(c(crossprod(O)))/4)/5
}
cyl <- cylinder3d(ctr, radius, sides = 30)
shade3d(cyl, color = "yellow")
}
# stereographic subdivision ####
midpoint4 <- function(A, B, r) {
xmid <- (A[1] + B[1]) / 2
ymid <- (A[2] + B[2]) / 2
zmid <- (A[3] + B[3]) / 2
tmid <- (A[4] + B[4]) / 2
mids <- c(xmid, ymid, zmid, tmid)
lg <- sqrt(c(crossprod(mids))) / r
mids / lg
}
subdiv0 <- function(A, B, C, r) {
mAB <- midpoint4(A, B, r)
mAC <- midpoint4(A, C, r)
mBC <- midpoint4(B, C, r)
trgl1 <- rbind(A, mAB, mAC)
trgl2 <- rbind(B, mBC, mAB)
trgl3 <- rbind(C, mAC, mBC)
trgl4 <- rbind(mAB, mAC, mBC)
list(trgl1, trgl2, trgl3, trgl4)
}
subdiv <- function(A, B, C, r, depth) {
if(depth == 1) {
out <- subdiv0(A, B, C, r)
} else {
triangles <- subdiv(A, B, C, r, depth-1)
out <- vector("list", 4^(depth-1))
for(i in 1:(4^(depth-1))) {
trgl <- triangles[[i]]
trgls <- subdiv0(trgl[1, ], trgl[2, ], trgl[3,], r)
out[[4*(i-1)+1]] = trgls[[1]]
out[[4*(i-1)+2]] = trgls[[2]]
out[[4*(i-1)+3]] = trgls[[3]]
out[[4*(i-1)+4]] = trgls[[4]]
}
}
out
}
#############################
#### draw the polychoron ####
#############################
theta <- pi/2
phi <- 0
r2 <- 3 + (1+sqrt(2))^2
r <- sqrt(r2)
depth <- 4
Draw <- function(xi) {
vs <- Vertices3(theta, phi, xi, r2)
stereoTriangles <- vector("list", 64)
for(i in 1:64) {
triangles4 <- subdiv(
vertices[faces[i, 1], ],
vertices[faces[i, 2], ],
vertices[faces[i, 3], ],
r, depth
)
stereoTriangles[[i]] <- vector("list", length(triangles4))
for(j in 1:length(triangles4)) {
trgl4 <- triangles4[[j]]
stereoTriangles[[i]][[j]] <-
rbind(
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[1, ]), r2
),
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[2, ]), r2
),
stereographicProjection(
rotate4d(theta, phi, xi, trgl4[3, ]), r2
)
)
}
}
# draw triangles
for(i in 1:64) {
trgl <- stereoTriangles[[i]]
vrtcs <- t(do.call(rbind, trgl))
indices <- matrix(1L:(3L*length(trgl)), nrow = 3L)
mesh <- tmesh3d(
vertices = vrtcs,
indices = indices
)
mesh <- Rvcg::vcgClean(mesh, sel = c(0, 7), silent = TRUE)
shade3d(mesh, color = "red")
}
# draw edges
for(i in 1:96) {
Edge(vertices, edges[i, 1], edges[i, 2], theta, phi, xi, r2)
}
# draw vertices
for(i in 1:64) {
spheres3d(
rbind(vs[i,]),
radius = log1p(sqrt(c(crossprod(vs[i, ])))/4)/4,
color = "yellow"
)
}
}
# plot ####
xi_ <- seq(0, 2*pi, length.out = 121)[-1]
open3d(windowRect = 50 + c(0, 0, 512, 512))
bg3d(
sphere = FALSE, texture = "SpaceBackground.png", col = "white"
)
view3d(0, 0, zoom = 0.9)
for(i in seq_along(xi_)) {
Draw(xi_[i])
snapshot3d(sprintf("zzpic%03d.png", i), webshot = FALSE)
clear3d()
}
pngs <- list.files(".", pattern = "^zzpic", full.names = TRUE)
library(gifski)
gifski(pngs, "runcinatedTesseract2_R.gif",
width = 512, height = 512, delay = 1/9)
file.remove(pngs)
@stla
Copy link
Author

stla commented Aug 5, 2023

runcinatedTesseract_tetrahedra

@stla
Copy link
Author

stla commented Aug 5, 2023

runcinatedTesseract2_R

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment