Solid Möbius strip to torus

SolidMobiusStrip_to_Torus.md

SolidMobiusStrip_to_Torus.R

library(misc3d)
library(rgl)
library(spray)

f <- function(x, y, z, a, b){
  ((x*x+y*y+1)*(a*x*x+b*y*y)+z*z*(b*x*x+a*y*y)-2*(a-b)*x*y*z-a*b*(x*x+y*y))^2 - 
    4*(x*x+y*y)*(a*x*x+b*y*y-x*y*z*(a-b))^2
}

gradient <- function(xyz){
  cbind(dfx(xyz), dfy(xyz), dfz(xyz))
}

# run the marching cubes algorithm ####
nx <- 220; ny <- 220; nz <- 220
x <- seq(-1.7, 1.7, length.out = nx) 
y <- seq(-1.7, 1.7, length.out = ny)
z <- seq(-0.7, 0.7, length.out = nz) 
G <- expand.grid(x = x, y = y, z = z)

b_ <- seq(0.05, 0.25, length.out = 50)

for(i in seq_along(b_)){
  cat("i: ", i, "\n")
  b <- b_[i]
  a <- 0.5 - b
  voxel <- array(with(G, f(x, y, z, a, b)), c(nx, ny, nz))
  surf <- computeContour3d(voxel, maxvol = max(voxel), level = 0, 
                           x = x, y = y, z = z)
  # build the rgl mesh ####
  P <- f(lone(1,3), lone(2,3), lone(3,3), a, b)
  dfx <- as.function(deriv(P, 1))
  dfy <- as.function(deriv(P, 2))
  dfz <- as.function(deriv(P, 3))
  mesh0 <- misc3d:::t2ve(makeTriangles(surf))
  mesh <- tmesh3d(
    vertices = mesh0$vb,
    indices = mesh0$ib,
    normals = -gradient(t(mesh0$vb))
  )
  # plot ####
  open3d(windowRect = c(50, 50, 562, 562))
  bg3d(rgb(54, 57, 64, maxColorValue = 255))
  view3d(0, -40, zoom = 0.8)
  shade3d(mesh, color="darkmagenta")
  rgl.snapshot(sprintf("zzpic%03d.png", i))
  close3d()
}

command <-
  "magick convert -dispose previous -delay 1x6 -duplicate 1,-2-1 zzpic*.png SolidMobiusStrip_to_Torus.gif"
system(command)