Soddy hexlet and Villarceau circles in R with rgl

SoddyHexletWithVillarceauCircles_rgl.R

# circumsphere ####
circumsphere <- function(p1, p2, p3, p4) {
  x1 <- p1[1]; y1 <- p1[2]; z1 <- p1[3]
  x2 <- p2[1]; y2 <- p2[2]; z2 <- p2[3]
  x3 <- p3[1]; y3 <- p3[2]; z3 <- p3[3]
  x4 <- p4[1]; y4 <- p4[2]; z4 <- p4[3]
  a <- det(cbind(rbind(p1, p2, p3, p4), 1))
  q1 <- c(crossprod(p1))
  q2 <- c(crossprod(p2))
  q3 <- c(crossprod(p3))
  q4 <- c(crossprod(p4))
  q <- c(q1, q2, q3, q4)
  x <- c(x1, x2, x3, x4)
  y <- c(y1, y2, y3, y4)
  z <- c(z1, z2, z3, z4)
  Dx <- det(cbind(q, y, z, 1))
  Dy <- -det(cbind(q, x, z, 1))
  Dz <- det(cbind(q, x, y, 1))
  c <- det(cbind(q, x, y, z))
  r <- 0.5*sqrt(Dx*Dx + Dy*Dy + Dz*Dz - 4*a*c) / abs(a) # ou distance(center,p1)
  list(center = 0.5*c(Dx, Dy, Dz)/a, radius=r)
}

# pretty rgl spheres ####
prettySpheres3d <- function(centers, radii, depth = 4, ...) {
  sphere <- subdivision3d(icosahedron3d(), depth=depth)
  sphere$vb[4, ] <- apply(sphere$vb[1:3, ], 2, function(x) sqrt(sum(x*x)))
  sphere$normals <- sphere$vb
  if(!is.matrix(centers)) {
    centers <- t(centers)
  }
  if(length(radii) == 1){
    radii <- rep(radii, nrow(centers))
  }
  for(i in 1:length(radii)) {
    shade3d(
      translate3d(
        scale3d(sphere, radii[i], radii[i], radii[i]), 
        centers[i, 1], centers[i, 2], centers[i, 3]), ...)
  }
  return(invisible())  
}

# torus passing by three points ####
plane3pts <- function(p1 ,p2, p3) { # plane passing by points p1, p2, p3 
  xcoef = (p1[2]-p2[2])*(p2[3]-p3[3])-(p1[3]-p2[3])*(p2[2]-p3[2]);
  ycoef = (p1[3]-p2[3])*(p2[1]-p3[1])-(p1[1]-p2[1])*(p2[3]-p3[3]);
  zcoef = (p1[1]-p2[1])*(p2[2]-p3[2])-(p1[2]-p2[2])*(p2[1]-p3[1]);
  offset = p1[1]*xcoef + p1[2]*ycoef + p1[3]*zcoef; 
  c(xcoef, ycoef, zcoef, offset)
}
circleCenterAndRadius <- function(p1, p2, p3) {
  # center, radius and normal of the circle passing by three points 
  p12 <- (p1+p2)/2
  p23 <- (p2+p3)/2
  v12 <- p2-p1
  v23 <- p3-p2
  plane1 <- plane3pts(p1, p2, p3)
  A <- rbind(plane1[1:3], v12, v23)
  b <- c(plane1[4], sum(p12*v12), sum(p23*v23))
  center <- c(solve(A) %*% b)
  r <- sqrt(c(crossprod(p1-center)))
  list(center = center, radius = r, normal = plane1[1:3])
}
transfoMatrix <- function(p1, p2, p3){ # transformation matrix 
  crn <- circleCenterAndRadius(p1, p2, p3)
  center <- crn$center
  radius <- crn$radius
  normal <- crn$normal
  if(normal[1]==0 && normal[2]==0){
    m <- rbind(cbind(diag(3), center), c(0, 0, 0, 1))
  } else {
    measure <- sqrt(normal[1]^2 + normal[2]^2 + normal[3]^2)
    normal <- normal/measure
    s <- sqrt(normal[1]^2 + normal[2]^2)
    u <- c(normal[2]/s, -normal[1]/s, 0)
    v <- geometry::extprod3d(normal, u)
    m <- rbind(cbind(u, v, normal, center), c(0, 0, 0, 1))
  }
  list(matrix = t(m), radius=radius)
}
createTorusMesh <- function(R, r, S, s, arc = 2*pi){
  vertices <- indices <- NULL
  for(j in 0:s) {
    v <- j / s * 2*pi
    for(i in 0:S) {
      u <- i / S * arc
      vertex <- c(
        (R + r*cos(v)) * cos(u),
        (R + r*cos(v)) * sin(u),
        r * sin(v)
      )
      vertices <- cbind(vertices, vertex)
    }
  }
  for(j in 1:s) {
    for(i in 1:S) {
      a <- (S + 1) * j + i 
      b <- (S + 1) * (j - 1) + i 
      c <- (S + 1) * (j - 1) + i + 1
      d <- (S + 1) * j + i + 1
      indices <- cbind(indices, c(a, b, d), c(b, c, d))
    }
  }
  tmesh3d(rbind(vertices, 1), indices)
}
torusMesh3points <- function(p1, p2, p3, r, S = 64, s=  32, arc = 2*pi) {
  mr <- transfoMatrix(p1, p2, p3)
  tmesh0 <- createTorusMesh(mr$radius, r, S, s, arc)
  transform3d(tmesh0, mr$matrix)
}

# parametrization of Villarceau circles ####
villarceau <- function(mu, a, c, theta, psi, epsilon) {
  b <- sqrt(a*a-c*c)
  bb <- b*sqrt(mu*mu-c*c)
  bb2 <- b*b*(mu*mu-c*c)
  denb1 <- c*(a*c-mu*c+c*c-a*mu-bb)
  b1 <- (a*mu*(c-mu)*(a+c)-bb2+c*c+bb*(c*(a-mu+c)-2*a*mu))/denb1
  denb2 <- c*(a*c-mu*c-c*c+a*mu+bb)
  b2 <- (a*mu*(c+mu)*(a-c)+bb2-c*c+bb*(c*(a-mu-c)+2*a*mu))/denb2
  omegaT <- (b1+b2)/2
  d <- (a-c)*(mu-c)+bb
  r <- c*c*(mu-c)/((a+c)*(mu-c)+bb)/d
  R <- c*c*(a-c)/((a-c)*(mu+c)+bb)/d
  omega2 <- (a*mu + bb)/c
  sign <- ifelse(epsilon > 0, 1, -1)
  f1 <- -sqrt(R*R-r*r)*sin(theta)
  f2 <- sign*(r+R*cos(theta))
  x1 <- f1*cos(psi) + f2*sin(psi) + omegaT
  y1 <- f1*sin(psi) - f2*cos(psi)
  z1 <- r*sin(theta)
  den <- (x1-omega2)^2 + y1*y1 + z1*z1
  c(omega2 + (x1-omega2)/den, y1/den, z1/den )
}

# Villarceau circle as mesh ####
vcircle <- function(mu, a, c, r, psi, sign, shift) {
  p1 <- villarceau(mu, a, c, 0, psi, sign) + shift
  p2 <- villarceau(mu, a, c, 2, psi, sign) + shift
  p3 <- villarceau(mu, a, c, 4, psi, sign) + shift
  torusMesh3points(p1, p2, p3, r)
}

# Soddy hexlet and its cyclide ####
inversion <- function(phi, point, radius, center) {
  omega <- c(radius/phi, 0, 0)
  k <- radius^2 * (1/phi^2 - 1)
  v <- point - omega - center
  omega + center - k/c(crossprod(v))*v
}
oneSphere <- function(phi, center, radius, beta) {
  pt <- center + c(radius*cos(beta)*2/3, radius*sin(beta)*2/3, 0)
  thRadius <- radius/3
  p1 <- pt + c(thRadius, 0, 0)
  p2 <- pt + c(0, thRadius, 0)
  p3 <- pt + c(-thRadius, 0, 0)
  p4 <- pt + c(0, 0, thRadius)
  q1 <- inversion(phi, p1, radius, center)
  q2 <- inversion(phi, p2, radius, center)
  q3 <- inversion(phi, p3, radius, center)
  q4 <- inversion(phi, p4, radius, center)
  cs <- circumsphere(q1, q2, q3, q4)
  list(center = cs$center - c(2*radius/phi, 0, 0), radius = cs$radius)
}
cyclide <- function(phi, center, radius) {
  thRadius <- radius/3
  p1 <- c(thRadius, 0, 0) 
  p2 <- c(0, thRadius, 0)
  p3 <- c(-thRadius, 0, 0)
  p4 <- c(0, 0, thRadius)
  q1 <- inversion(phi, p1, radius, c(0, 0, 0))
  q2 <- inversion(phi, p2, radius, c(0, 0, 0))
  q3 <- inversion(phi, p3, radius, c(0, 0, 0))
  q4 <- inversion(phi, p4, radius, c(0, 0, 0))
  cs <- circumsphere(q1, q2, q3, q4)
  O2 <- cs$center
  list(mu = (radius - cs$radius) / 2, 
       a = (radius + cs$radius) / 2,
       c = -radius/phi + O2[1]/2,
       shift = center + (O2 - c(2/phi*radius, 0, 0)) / 2)
}
hexlet <- function(phi, cr, frame = 0, clockwise = TRUE){
  shift <- frame * ifelse(clockwise, -pi/90, pi/90)
  sapply(
    c(0, pi/3, 2*pi/3,pi, 4*pi/3, 5*pi/3) + shift,  
    function(beta) oneSphere(phi, cr$center, cr$radius, beta)
  )
}


# rgl ####
phi <- 0.25
library(rgl)
hx <- hexlet(phi, list(center = c(0, 0, 0), radius = 1))
open3d(windowRect = c(50, 50, 562, 562))
view3d(0, 0)
par3d(zoom = 0.8)
centers <- t(apply(hx, 2, `[[`, 1)); radii <- apply(hx, 2, `[[`, 2)
prettySpheres3d(centers, radii, color = "black")
cycld <- cyclide(phi, c(0, 0, 0), 1)
n <- 15
for(i in 1:n) {
  psi <- i*2*pi/n
  vmesh <- vcircle(cycld$mu, cycld$a, cycld$c, 0.01, psi, 1, cycld$shift)
  shade3d(vmesh, color = "red")
  vmesh <- vcircle(cycld$mu, cycld$a, cycld$c, 0.01, psi, -1, cycld$shift)
  shade3d(vmesh, color = "red")
}
snapshot3d("SoddyHexletWithVillarceauCircle.png", webshot = FALSE)

# animation ####
library(rgl)
open3d(windowRect = c(50, 50, 562, 562))
view3d(0, 0)
par3d(zoom = 0.8)
phi <- 0.25
cycld <- cyclide(phi, c(0, 0, 0), 1)
n <- 15
for(frame in 1:30) {
  for(i in 1:n) {
    psi <- i*2*pi/n
    vmesh <- vcircle(cycld$mu, cycld$a, cycld$c, 0.01, psi, 1, cycld$shift)
    shade3d(vmesh, color = "red")
    vmesh <- vcircle(cycld$mu, cycld$a, cycld$c, 0.01, psi, -1, cycld$shift)
    shade3d(vmesh, color = "red")
  }
  hx <- hexlet(phi, list(center = c(0, 0, 0), radius = 1), frame = frame)
  quads3d(rbind(c(1, 1, -1),
                c(1, -1, -1),
                c(-1, -1, -1),
                c(-1, 1, -1)), 
          color = "white")
  prettySpheres3d(
    t(apply(hx, 2, `[[`, 1)), apply(hx, 2, `[[`, 2), color = "black"
  )
  snapshot3d(sprintf("SoddyVillarceau%02d.png", frame), webshot = FALSE)
  clear3d()
}

command <- 
  "magick convert -dispose previous SoddyVillarceau*.png -loop 0 -delay 5 -alpha on SoddyVillarceau.gif"
system(command)

pngs <- list.files(pattern = "SoddyVillarceau.*png")
file.remove(pngs)