Some implicit surfaces in R, with misc3d and rgl

.gitignore

.Rproj.user
.Rhistory
.RData
.Ruserdata
*.Rproj

BarthSextic1.R

phi <- (1+sqrt(5))/2
f0 <- function(x,y,z){
  4*(phi^2*x^2-y^2)*(phi^2*y^2-z^2)*(phi^2*z^2-x^2) - 
    (1+2*phi)*(x^2+y^2+z^2-1)^2*1
}
f1 <- function(xyz){
  x <- xyz[1]; y <- xyz[2]; z <- xyz[3]
  f0(x,y,z)
}

library(numDeriv)
gradient <- function(xyz){
  grad(f1, xyz)
}

nx <- 100; ny <- 100; nz <- 100
x <- seq(-2, 2, length=nx) 
y <- seq(-2, 2, length=ny)
z <- seq(-2, 2, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f0(x,y,z)), c(nx,ny,nz))

library(misc3d)
surf <- computeContour3d(voxel, level=0, x=x, y=y, z=z)

library(rgl)
mesh <- tmesh3d(vertices = t(surf), 
                indices = matrix(1:nrow(surf), nrow=3), 
                homogeneous = FALSE, 
                normals = -t(apply(surf, 1, gradient)))

open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color=rgb(1,0,1))

BarthSextic2.R

phi <- (1+sqrt(5))/2
f0 <- function(x,y,z){
  4*(phi^2*x^2-y^2)*(phi^2*y^2-z^2)*(phi^2*z^2-x^2) - 
    (1+2*phi)*(x^2+y^2+z^2-1)^2*1
}
gradient <- function(xyz){
  x <- xyz[,1]; y <- xyz[,2]; z <- xyz[,3]  
  x2 <- x*x; y2 <- y*y; z2 <- z*z; phi2 <- phi*phi
  cbind(
    8*x*phi2*(z2*phi2-x2)*(y2*phi2-z2) - 8*x*(x2*phi2-y2)*(y2*phi2-z2) -
      4*x*(2*phi+1)*(x2+y2+z2-1),
    8*y*phi2*(x2*phi2-y2)*(z2*phi2-x2) - 8*y*(z2*phi2-x2)*(y2*phi2-z2) -
      4*y*(2*phi+1)*(x2+y2+z2-1),
    8*z*phi2*(x2*phi2-y2)*(y2*phi2-z2) - 8*z*(x2*phi2-y2)*(z2*phi2-x2) -
      4*z*(2*phi+1)*(x2+y2+z2-1)
  )
}

nx <- 120; ny <- 120; nz <- 120
x <- seq(-1.8, 1.8, length=nx) 
y <- seq(-1.8, 1.8, length=ny)
z <- seq(-1.8, 1.8, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f0(x,y,z)), c(nx,ny,nz))
mask <- array(with(g, x^2+y^2+z^2<3), c(nx,ny,nz))
voxel[!mask] <- NaN # <=> using mask=mask in computeContour3d
library(misc3d)
surf <- computeContour3d(voxel, maxvol=max(voxel[mask]), level=0, x=x, y=y, z=z)

library(rgl)
mesh <- tmesh3d(vertices = t(surf),
                indices = matrix(1:nrow(surf), nrow=3),
                homogeneous = FALSE,
                normals = -gradient(surf))

open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color=rgb(1,0,1))

C8surface.R

f <- function(x,y,z){
  64*x^8 - 128*x^6 + 80*x^4 - 16*x^2 + 2 + 64*y^8 - 128*y^6 + 80*y^4 - 16*y^2 + 
    64*z^8 - 128*z^6 + 80*z^4 - 16*z^2
}
gradient <- function(xyz){
  x <- xyz[1]; y <- xyz[2]; z <- xyz[3]
  c(
    8*64*x^7 - 6*128*x^5 + 4*80*x^3 - 2*16*x,
    8*64*y^7 - 6*128*y^5 + 4*80*y^3 - 2*16*y,
    8*64*z^7 - 6*128*z^5 + 4*80*z^3 - 2*16*z
  )
}

nx <- 100; ny <- 100; nz <- 100
x <- seq(-1.1, 1.1, length=nx) 
y <- seq(-1.1, 1.1, length=ny)
z <- seq(-1.1, 1.1, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f(x,y,z)), c(nx,ny,nz))

library(misc3d)
surf <- computeContour3d(voxel, level=0, x=x, y=y, z=z)

library(rgl)
mesh <- tmesh3d(vertices = t(surf), 
                indices = matrix(1:nrow(surf), nrow=3), 
                homogeneous = FALSE, 
                normals = -t(apply(surf, 1, gradient)))

open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color=rgb(1,0,1))

Decocube.R

f <- function(x,y,z,a){
  x2 = x*x
  y2 = y*y
  z2 = z*z
  a2 = a*a
  x2y2a2 = x2+y2-a2
  y2z2a2 = y2+z2-a2
  z2x2a2 = z2+x2-a2
  (x2y2a2*x2y2a2 + (z2-1)*(z2-1)) *
    (y2z2a2*y2z2a2 + (x2-1)*(x2-1)) *
    (z2x2a2*z2x2a2 + (y2-1)*(y2-1))
}

gradient <- function(xyz,a){
  x <- xyz[1]; y <- xyz[2]; z <- xyz[3]
  x2 = x*x
  y2 = y*y
  z2 = z*z
  a2 = a*a
  x2y2a2 = x2+y2-a2
  y2z2a2 = y2+z2-a2
  z2x2a2 = z2+x2-a2
  c(
    4*x*(x2-1)*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)) + 
      4*x*x2y2a2*(y2z2a2*y2z2a2 + (x2-1)*(x2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)) + 
      4*x*z2x2a2*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (y2z2a2*y2z2a2 + (x2-1)*(x2-1)),
    4*y*(y2-1)*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (y2z2a2*y2z2a2 + (x2-1)*(x2-1)) + 
      4*y*x2y2a2*(y2z2a2*y2z2a2 + (x2-1)*(x2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)) + 
      4*y*y2z2a2*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)),
    4*z*(z2-1)*(y2z2a2*y2z2a2 + (x2-1)*(x2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)) + 
      4*z*y2z2a2*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (z2x2a2*z2x2a2 + (y2-1)*(y2-1)) + 
      4*z*z2x2a2*(x2y2a2*x2y2a2 + (z2-1)*(z2-1)) * 
      (y2z2a2*y2z2a2 + (x2-1)*(x2-1))
  )
}

a = 0.95

nx <- 150; ny <- 150; nz <- 150
x <- seq(-1.3, 1.3, length=nx) 
y <- seq(-1.3, 1.3, length=ny)
z <- seq(-1.3, 1.3, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f(x,y,z,a)), c(nx,ny,nz))

library(misc3d)
surf <- computeContour3d(voxel, level=0.01, x=x, y=y, z=z)

library(rgl)
mesh <- tmesh3d(vertices = t(surf), 
                indices = matrix(1:nrow(surf), nrow=3), 
                homogeneous = FALSE, 
                normals = -t(apply(surf, 1, function(xyz){
                  gradient(xyz,a)
                })))
open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color="whitesmoke")

Enzensberger.R

# Enzensberger star #

library(misc3d)
library(rgl)

f <- function(x,y,z){
  400*(x^2*y^2+ y^2*z^2+ x^2*z^2)- (1-x^2-y^2-z^2)^3 
}

gradient <- function(xyz){
  x <- xyz[,1]; y <- xyz[,2]; z <- xyz[,3]
  cbind(
    6*x*(-x^2 - y^2 - z^2 + 1)^2 + 800*x*(y^2 + z^2),
    6*y*(-x^2 - y^2 - z^2 + 1)^2 + 800*y*(x^2 + z^2),
    6*z*(-x^2 - y^2 - z^2 + 1)^2 + 800*z*(x^2 + y^2)
  )
}

n <- 150
x <- y <- z <- seq(-0.9, 0.9, len=n)
voxel <- array(with(expand.grid(x=x, y=y, z=z), f(x,y,z)), dim=c(n,n,n))
surf <- computeContour3d(voxel, max(voxel), 0, x=x, y=y, z=z)

mesh <- tmesh3d(vertices = t(surf), 
                indices = matrix(1:nrow(surf), nrow=3), 
                homogeneous = FALSE, 
                normals = -gradient(surf))

open3d(windowRect = c(50,50,550,550))
bg3d(rgb(54,57,64, maxColorValue = 255))
shade3d(mesh, color = "orange")

# animation ####
movie3d(spin3d(axis = c(0, 0, 1), rpm = 60),
        duration = 1, fps = 60,
        movie = "Enzensberger", dir = ".",
        convert = "convert -dispose previous -loop 0 -delay 1x%d %s*.png %s.%s",
        startTime = 1/60)

ICN5D_cage.R

library(rgl)
library(misc3d)
library(Rvcg)

a <- 0

# f <- function(x, y, z, a, t){
#   (((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)^2-4*(9*z^2+16*(x*sin(t)+a*cos(t))^2))^2*(((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)^2+296*((x*cos(t)-a*sin(t))^2+y^2)-4*(9*z^2+16*(x*sin(t)+a*cos(t))^2)) -16*((x*cos(t)-a*sin(t))^2+y^2)*((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)^2*(37*((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)^2-1369*((x*cos(t)-a*sin(t))^2+y^2)-7*(225*z^2+448*(x*sin(t)+a*cos(t))^2)) -16*sqrt(3)/9*((x*cos(t)-a*sin(t))^3-3*(x*cos(t)-a*sin(t))*y^2)*(110*((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)^3 -148*((x*cos(t)-a*sin(t))^2+y^2+z^2+(x*sin(t)+a*cos(t))^2+145/3)*(110*(x*cos(t)-a*sin(t))^2+110*y^2-297*z^2+480*(x*sin(t)+a*cos(t))^2)) -64*((x*cos(t)-a*sin(t))^2+y^2)*(3*(729*z^4+4096*(x*sin(t)+a*cos(t))^4)+168*((x*cos(t)-a*sin(t))^2+y^2)*(15*z^2-22*(x*sin(t)+a*cos(t))^2)) +64*(12100/27*((x*cos(t)-a*sin(t))^3-3*(x*cos(t)-a*sin(t))*y^2)^2 -7056*(3*(x*cos(t)-a*sin(t))^2*y-y^3)^2) -592240896*z^2*(x*sin(t)+a*cos(t))^2
# }
# since a = 0
f <- function(x, y, z, t){
  (((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)^2-4*(9*z^2+16*(x*sin(t))^2))^2*(((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)^2+296*((x*cos(t))^2+y^2)-4*(9*z^2+16*(x*sin(t))^2)) -16*((x*cos(t))^2+y^2)*((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)^2*(37*((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)^2-1369*((x*cos(t))^2+y^2)-7*(225*z^2+448*(x*sin(t))^2)) -16*sqrt(3)/9*((x*cos(t))^3-3*(x*cos(t))*y^2)*(110*((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)^3 -148*((x*cos(t))^2+y^2+z^2+(x*sin(t))^2+145/3)*(110*(x*cos(t))^2+110*y^2-297*z^2+480*(x*sin(t))^2)) -64*((x*cos(t))^2+y^2)*(3*(729*z^4+4096*(x*sin(t))^4)+168*((x*cos(t))^2+y^2)*(15*z^2-22*(x*sin(t))^2)) +64*(12100/27*((x*cos(t))^3-3*(x*cos(t))*y^2)^2 -7056*(3*(x*cos(t))^2*y-y^3)^2) -592240896*z^2*(x*sin(t))^2
}

nx <- 250; ny <- 250; nz <- 250
x <- seq(-8, 9, length.out = nx)
y <- seq(-9, 9, length.out = ny)
z <- seq(-6, 6, length.out = nz)
G <- expand.grid(x = x, y = y, z = z)

userMatrix <- rbind(
  c( 0.69, -0.71, 0.12, 0),
  c( 0.00,  0.16, 0.99, 0),
  c(-0.72, -0.68, 0.11, 0),
  c( 0.00,  0.00, 0.00, 1)
)

t_ <- seq(0, pi/2, length.out = 50)

for(i in seq_along(t_)){
  # run the marching cubes algorithm ####
  voxel <- array(with(G, f(x, y, z, t_[i])), dim = c(nx, ny, nz))
  surface <- computeContour3d(
    voxel, maxvol = max(voxel), level = 0, x = x, y = y, z = z
  )
  # make the mesh
  mesh0 <- misc3d:::t2ve(makeTriangles(surface))
  mesh <- vcgUpdateNormals(tmesh3d(
    vertices = mesh0$vb,
    indices = mesh0$ib
  ))
  # plot
  open3d(windowRect = c(50, 50, 562, 562))
  bg3d(rgb(54, 57, 64, maxColorValue = 255))
  par3d(userMatrix = userMatrix, zoom = 0.75)
  shade3d(mesh, color = "turquoise4")
  rgl.snapshot(sprintf("zzzpic%03d.png", i))
  close3d()
}

# animation
command <-
  "magick convert -dispose previous -delay 1x8 -duplicate 1,-2-1 zzzpic*.png ICN5Dcage.gif"
system(command)

pngs <- list.files(".", pattern = "^zzzpic", full.names = TRUE)
file.remove(pngs)

ICN5D_eight.R

library(rgl)
library(misc3d)

a <- 2.4

f <- function(x, y, z, a){
  (sqrt((x^2 - a^2)^2/(x^2 + a^2) + y^2) - 3)^2 +
    (sqrt((x*a)^2/(x^2 + a^2) + z^2) - 1.5)^2
}

# run the marching cubes algorithm ####
nx <- 200; ny <- 200; nz <- 200
x <- seq(-5, 5, length.out = nx)
y <- seq(-4, 4, length.out = ny)
z <- seq(-2.5, 2.5, length.out = nz)
G <- expand.grid(x = x, y = y, z = z)
voxel <- array(with(G, f(x, y, z, a)), dim = c(nx, ny, nz))
surface <- computeContour3d(
  voxel, maxvol = max(voxel), level = 0.5, x = x, y = y, z = z
)

mesh0 <- misc3d:::t2ve(makeTriangles(surface))
mesh <- addNormals(tmesh3d(
  vertices = mesh0$vb,
  indices = mesh0$ib
))

open3d(windowRect = c(50, 50, 562, 562))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
view3d(30, -30, zoom= 0.75)
shade3d(mesh, color = "tomato")

# 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 = "echo \"%d %s %s %s\"",
  clean = FALSE,
  webshot = FALSE
)

pngs <- list.files(".", pattern = "^zzzpic", full.names = TRUE)
library(gifski)
gifski(pngs, "ICN5D_eight.gif",
       width = 512, height = 512, delay = 1/9)

file.remove(pngs)

ImplicitSurfaces.md

KohnNirenberg.R

library(misc3d)
library(spray)

f <- function(x, y, z, k){
  b <- x^2 + y^2
  z <- 3*z
  z + z^2*b + b^4 + 15/7*b*(x^6 - 15*x^4*y^2 + 15*x^2*y^4 - y^6) + 
    k*(x^2 + y^2 + z^2)^5
}

# f as a polynomial
P <- f(lone(1,3), lone(2,3), lone(3,3), 0.1)

# gradient of f
dfx <- as.function(deriv(P, 1))
dfy <- as.function(deriv(P, 2))
dfz <- as.function(deriv(P, 3))
gradient <- function(xyz){
  cbind(dfx(xyz), dfy(xyz), dfz(xyz))
}

# compute isosurface
n <- 150
x <- seq(-3.1, 3.1, length.out = n)
y <- seq(-3.5, 3.5, length.out = n)
z <- seq(-0.6, 0.5, length.out = n)
g <- expand.grid(x=x, y=y, z=z)
v <- array(with(g, f(x,y,z,k=0.1)), dim = c(n,n,n))
surf <- computeContour3d(v, max(v), 0, x=x, y=y, z=z, mask = NULL)

# make rgl mesh
library(rgl)
mesh <- tmesh3d(vertices = t(surf),
                indices = matrix(1:nrow(surf), nrow=3),
                homogeneous = FALSE,
                normals = -gradient(surf))

# draw
open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
view3d(0, -30, zoom = 0.75)
shade3d(mesh, color="gold")

# animation ####
movie3d(spin3d(axis = c(0, 0, 1), rpm = 60),
        duration = 1, fps = 60,
        movie = "KohnNirenberg00", dir = ".",
        convert = "magick convert -dispose previous -loop 0 -delay 1x%d %s*.png %s.%s",
        startTime = 1/60)

LabsCubic.R

library(misc3d)
library(rgl)

f <- function(ρ, θ, ϕ){
  x <- ρ*cos(θ)*sin(ϕ)
  y <- ρ*sin(θ)*sin(ϕ)
  z <- ρ*cos(ϕ)
  x^3 + 3*x^2 - 3*x*y^2 + 3*y^2 - 4 + z^3 + 3*z^2
}

nρ <- 200; nθ <- 200; nϕ <- 200
ρ <- seq(0, 6, length = nρ) 
θ <- seq(0, 2*pi, length = nθ)
ϕ <- seq(0, pi, length = nϕ) 
g <- expand.grid(ρ=ρ, θ=θ, ϕ=ϕ)
voxel <- array(with(g, f(ρ,θ,ϕ)), c(nρ,nθ,nϕ))

surf <- computeContour3d(voxel, maxvol = max(voxel), level = 0)
ρ <- (surf[,1]-1) * 6 / (nρ - 1)
θ <- (surf[,2]-1) * 2*pi / (nθ - 1)
ϕ <- (surf[,3]-1) * pi / (nϕ - 1)
surf <- cbind(ρ*cos(θ)*sin(ϕ), ρ*sin(θ)*sin(ϕ), ρ*cos(ϕ))

gradient <- function(xyz){
  x <- xyz[,1]; y <- xyz[,2]; z <- xyz[,3]
  cbind(
    3*x^2 + 6*x - 3*y^2,
    -6*x*y + 6*y,
    3*z^2 + 6*z
  )
}

mesh <- tmesh3d(vertices = t(surf),
                indices = matrix(1:nrow(surf), nrow=3),
                homogeneous = FALSE,
                normals = gradient(surf))

# make colors
distances <- sqrt(apply(surf, 1, crossprod))
dmin <- min(distances); dmax <- max(distances)
fcolor <- function(d, palette){
  clrRamp <- colorRamp(palette)
  RGB <- clrRamp((d-dmin)/(dmax-dmin))
  rgb(RGB[,1], RGB[,2], RGB[,3], maxColorValue = 255)
}
colors1 <- fcolor(distances, jcolors::jcolors_contin("pal10", reverse = TRUE)(256))
colors2 <- fcolor(distances, jcolors::jcolors_contin("pal2")(256))

# draw
open3d(windowRect = c(50, 50, 550, 550), zoom = 0.8)
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color = colors1, meshColor = "vertices", back = "culled")
shade3d(mesh, color = colors2, meshColor = "vertices", front = "culled")

# animation #### 
movie3d(spin3d(axis = c(0, 0, 1), rpm = 10),
        duration = 6, fps = 20,
        movie = "LabsCubic", dir = ".",
        convert = "gifski --fps %d %s*.png -o %s.%s",
        startTime = 1/120)

NordstrandTwenty.R

f0 <- function(x,y,z){
  25*(x^3*(y+z)+y^3*(x+z)+z^3*(x+y)) + 50*(x^2*y^2+x^2*z^2+y^2*z^2) - 
    125*(x^2*y*z+y^2*x*z+z^2*x*y) + 60*x*y*z - 4*(x*y+x*z+y*z)
}

gradient <- function(xyz){
  x <- xyz[,1]; y <- xyz[,2]; z <- xyz[,3]
  cbind(
    25*(3*x^2*(y+z) + y^3 + z^3) + 100*x*(y^2 + z^2) - 125*y*z*(2*x + y + z) + 
      60*y*z - 4*(y + z),
    25*(3*y^2*(x+z) + x^3 + z^3) + 100*y*(x^2 + z^2) - 125*x*z*(2*y + x + z) + 
      60*x*z - 4*(x + z),
    25*(3*z^2*(x+y) + x^3 + y^3) + 100*z*(x^2 + y^2) - 125*x*y*(2*z + x + y) + 
      60*x*y - 4*(x + y)
  )
}

nx <- 120; ny <- 120; nz <- 120
x <- seq(-1.8, 1.8, length=nx) 
y <- seq(-1.8, 1.8, length=ny)
z <- seq(-1.8, 1.8, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f0(x,y,z)), c(nx,ny,nz))
mask <- array(with(g, x^2+y^2+z^2<3), c(nx,ny,nz))
voxel[!mask] <- NaN 
library(misc3d)
surf <- computeContour3d(voxel, maxvol=max(voxel[mask]), level=0, x=x, y=y, z=z)

library(rgl)
mesh0 <- tmesh3d(vertices = t(surf), 
                 indices = matrix(1:nrow(surf), nrow=3), 
                 homogeneous = FALSE, 
                 normals = -gradient(surf))

# "flower" (twenty Nordstrand surfaces) ####
source("ConeMesh.R") # https://gist.github.com/stla/b3db4a593c9b030567bde701d45fa385

Reorient_Trans <- function(Axis1, Axis2){
  vX1 <- Axis1 / sqrt(c(crossprod(Axis1)))
  vX2 <- Axis2 / sqrt(c(crossprod(Axis2)))
  Y <- rgl:::xprod(vX1, vX2)
  vY <- Y / sqrt(c(crossprod(Y)))
  Z1 <- rgl:::xprod(vX1, vY)
  vZ1 <- Z1 / sqrt(c(crossprod(Z1)))
  Z2 <- rgl:::xprod(vX2, vY)
  vZ2 <- Z2 / sqrt(c(crossprod(Z2)))
  M1 <- cbind(vX1, vY, vZ1)
  M2 <- rbind(vX2, vY, vZ2)
  M <- M1 %*% M2
  rbind(cbind(M, c(0,0,0)), c(0,0,0,1))
}

# vertices ####
phi <- (1+sqrt(5))/2
a <- 1/sqrt(3)
b <- a/phi 
c <- a*phi
vertices <- 
  rbind(
    c( a,  a,  a), 
    c( a,  a, -a),
    c( a, -a,  a),
    c(-a, -a,  a),
    c(-a,  a, -a),
    c(-a,  a,  a),
    c( 0,  b, -c), 
    c( 0, -b, -c), 
    c( 0, -b,  c),
    c( c,  0, -b),
    c(-c,  0, -b),
    c(-c,  0,  b),
    c( b,  c,  0),
    c( b, -c,  0),
    c(-b, -c,  0),
    c(-b,  c,  0),
    c( 0,  b,  c),
    c( a, -a, -a),
    c( c,  0,  b),
    c(-a, -a, -a)
  )

# rgl ####
colors <- randomcoloR::distinctColorPalette(20)
mesh0$normals <- rbind(mesh0$normals, 1)
sc <- 0.175
mesh <- rotate3d(
  rotate3d(
    scale3d(mesh0, sc, sc, sc), 
    pi/4, 0, 1, 0),
  -35*pi/180, 1, 0, 0)

open3d(windowRect=c(50,50,550,550))
bg3d(rgb(54,57,64, maxColorValue = 255))
for(i in 1:20){
  v <- vertices[i,]
  M <- Reorient_Trans(c(0,0,-1), v)
  obj <- transform3d(mesh, M)
  obj <- translate3d(obj, v[1], v[2], v[3])
  shade3d(obj, color=colors[i])
  shade3d(ConeMesh(c(0,0,0), 0, v, 0.08, 6, 40, color="darkolivegreen"))
}

Orthocircle.R

f <- function(x,y,z,a,b){
  x2 <- x*x
  y2 <- y*y
  z2 <- z*z
  xy2 <- x2+y2-1
  yz2 <- y2+z2-1
  zx2 <- z2+x2-1
  (xy2*xy2 + z2) * (yz2*yz2 + x2) * (zx2*zx2 + y2) - a*a*(1 + b*(x2 + y2 + z2))
}

gradient <- function(xyz,a,b){
  x <- xyz[1]; y <- xyz[2]; z <- xyz[3]
  x2 <- x*x
  y2 <- y*y
  z2 <- z*z
  xy2 <- x2+y2-1
  yz2 <- y2+z2-1
  zx2 <- z2+x2-1
  c(
    -2*a*a*b*x + 2*x*(xy2*xy2+z2)*(zx2*zx2+y2) + 
      4*x*zx2*(xy2*xy2+z2)*(yz2*yz2+x2) + 
      4*x*xy2*(zx2*zx2+y2)*(yz2*yz2+x2),
    -2*a*a*b*y + 2*y*(xy2*xy2+z2)*(yz2*yz2+x2) + 
      4*y*yz2*(xy2*xy2+z2)*(zx2*zx2+y2) + 
      4*y*xy2*(zx2*zx2+y2)*(yz2*yz2+x2),
    -2*a*a*b*z + 2*z*(zx2*zx2+y2)*(yz2*yz2+x2) + 
      4*z*yz2*(xy2*xy2+z2)*(zx2*zx2+y2) + 
      4*z*zx2*(xy2*xy2+z2)*(yz2*yz2+x2)
  )
}

a = 0.075; b = 3

nx <- 100; ny <- 100; nz <- 100
x <- seq(-1.3, 1.3, length=nx) 
y <- seq(-1.3, 1.3, length=ny)
z <- seq(-1.3, 1.3, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f(x,y,z,a,b)), c(nx,ny,nz))

library(misc3d)
surf <- computeContour3d(voxel, level=0, x=x, y=y, z=z)

library(rgl)
mesh <- tmesh3d(vertices = t(surf), 
                indices = matrix(1:nrow(surf), nrow=3), 
                homogeneous = FALSE, 
                normals = -t(apply(surf, 1, function(xyz){
                  gradient(xyz,a,b)
                })))

open3d(windowRect = c(50, 50, 550, 550))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
shade3d(mesh, color="whitesmoke")

SoccerStar.R

library(rgl)
library(rmarchingcubes)

# soccer star
# isosurface f=0
a1 <- a2 <- a3 <- a4 <- -100
f <- function(x, y, z){
  u <- x*x + y*y + z*z
  v <- -z * (2*x+z) *
    (x^4 - x*x*z*z + z^4 + 2*(x^3*z-x*z^3) + 5*(y^4-y*y*z*z) + 
       10*(x*y*y*z-x*x*y*y))
  w <- (4*x*x + z*z - 6*x*z) * 
    (z^4 - 2*z^3*x - x*x*z*z + 2*z*x^3 + x^4 - 25*y*y*z*z - 30*x*y*y*z - 
       10*x*x*y*y + 5*y^4) *
    (z^4 + 8*z^3*x + 14*x*x*z*z - 8*z*x^3 + x^4 - 10*y*y*z*z - 
       10*x*x*y*y + 5*y^4)
  1 + ((128565+115200*sqrt(5))/1295029 * a3 + 
         (49231296000*sqrt(5)-93078919125)/15386239549 * a4 - 
         a1 - 3*a2 - 3) * u +
    ((-230400*sqrt(5) - 257130)/1295029 * a3 + 
       (238926989250-126373248000*sqrt(5))/15386239549 * a4 + 3*a1 + 
       8*a2 + 3) * u * u + 
    ((115200*sqrt(5)+128565)/1295029 * a3 + 
       (91097280000*sqrt(5)-172232645625)/15386239549 * a4 - 3*a1 - 
       6*a2 - 1) * u*u*u + 
    (a3 + (121075-51200*sqrt(5))/11881 * a4) * v + 
    ((102400*sqrt(5)-242150)/11881 - 2*a3) * u * v + 
    a1 * u^4 + a2 * u^5 + a3 * u*u*v + a4*w
}

# make the isosurface
nx <- 400L; ny <- 300L; nz <- 300L
x <- seq(-1.5, 1.5, length.out = nx) 
y <- seq(-1.5, 1.5, length.out = ny)
z <- seq(-1.5, 1.5, length.out = nz) 
Grid <- expand.grid(X = x, Y = y, Z = z)
voxel <- array(with(Grid, f(X, Y, Z)), dim = c(nx, ny, nz))
cont <- contour3d(voxel, level = 0, x = x, y = y, z = z)

# plot
rmesh <- tmesh3d(
  vertices = t(cont[["vertices"]]),
  indices  = t(cont[["triangles"]]),
  normals  = cont[["normals"]],
  homogeneous = FALSE
)

open3d(windowRect = 50 + c(0, 0, 512, 512))
bg3d("#363940")
view3d(0, 0, zoom = 0.7)
shade3d(rmesh, color = "hotpink")

# 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, "SoccerStar.gif",
       width = 512, height = 512, 
       delay = 1/10)

file.remove(pngs)

decocube_metamorphosis.gif

decocube_metamorphosis.R

doubleTrumpet.gif

gluingTori.gif

gluingTori.R

icosahedralStar.R

implicitizationSatellite_k3.gif

implicitizationSatellite_k3.R

SoccerStar.gif

SolidMobiusStrip.R

SpaceBackground.png

SpiderCage.R

Tiger.R

Togliatti.R

Togliatti2.R

TogliattiTwenty.R

Valentine.gif

Valentine.R

Valentine2.gif

Valentine2.R

vanStraten.R

Icosahedral star