Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Recall that a Hopf torus is a two-dimensional object in the 4D space defined by a profile curve: a closed curve on the unit sphere. When mapping it to the 3D space with the stereographic projection, we can obtain a beautiful object, or an ugly one, depending on the choice of the profile curve.

Here, we will see the Hopf torus defined by this profile curve: This is the intersection of the unit sphere with a cubical cone, the isosurface of equation $$x^2 + y^2 + z^2 = 0$$. First, I obtained a mesh of this surface thanks to the rmarchingcubes package:

library(rgl)
library(rmarchingcubes)

# cubical cone ####
f <- function(x, y, z) {
x^3 + y^3 + z^3
}
x_ <- y_ <- z_ <- seq(-1.05, 1.05, length.out = 150L)
Grid <- expand.grid(X = x_, Y = y_, Z = z_)
voxel <- array(with(Grid, f(X, Y, Z)), dim = c(150L, 150L, 150L))
surf <- contour3d(voxel, level = 0, x_, y_, z_)
coneMesh <- tmesh3d(
vertices = t(surf[["vertices"]]),
indices  = t(surf[["triangles"]]),
normals  = surf[["normals"]]
)

Then I obtained the intersection with the unit sphere thanks to the clipMesh3d and getBoundary3d functions of the rgl package:

# intersection with unit sphere ####
sphereMesh <- Rvcg::vcgSphere(subdivision = 5L)
mesh <- clipMesh3d(
sphereMesh, fn = f, minVertices = 20000L
)
boundary <-
getBoundary3d(mesh, sorted = TRUE, color = "black", lwd = 2)
# plot ####
open3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.9)
shade3d(coneMesh, color = "red", alpha = 0.5, polygon_offset = 1)
shade3d(sphereMesh, color = "blue", alpha = 0.5, polygon_offset = 1)
shade3d(boundary)

The curve has not the required orientation of a nice profile curve. Its axis of symmetry is directed by $$(1,1,1)$$, and we need $$(0,0,1)$$ instead. So one has to rotate the curve. To do so, I use an exported function from the cgalMeshes package, namely quaternionFromTo. It will return a unit quaternion corresponding to the desired rotation. I already talked about this way to obtain a rotation sending a given vector to another given vector, here.

# get points at the intersection and rotate them ####
pts <- boundary[["vb"]][-4L, boundary[["is"]][1L, ]]
q <- cgalMeshes:::quaternionFromTo(c(1, 1, 1)/sqrt(3), c(0, 0, 1))
R <- onion::as.orthogonal(q)
gamma0 <- t(R %*% pts)[, c(3L, 2L, 1L)]

Now let’s introduce a function which creates a mesh of the Hopf torus defined by a discretized curve, such as our matrix of points gamma0. Again, I use an exported function from cgalMeshes, namely meshTopology, which returns the incidences between the vertices of the mesh.

# Hopf torus mesh from a discrete curve gamma ####
hMesh <- function(gamma, m) {
nu <- nrow(gamma)
uperiodic <- TRUE
u_ <- 1L:nu
vperiodic <- TRUE
nv <- as.integer(m)
v_ <- 1L:nv
R <- array(NA_real_, dim = c(3L, nv, nu))
for(k in 1L:nv) {
K <- k - 1L
cosphi <- cospi(2*K/m)
sinphi <- sinpi(2*K/m)
for(j in 1L:nu) {
p1 <- gamma[j, 1L]
p2 <- gamma[j, 2L]
p3 <- gamma[j, 3L]
yden <- sqrt(2 * (1 + p1))
y1 <- (1 + p1) / yden
y2 <- p2 / yden
y3 <- p3 / yden
x1 <- cosphi * y3 + sinphi * y2
x2 <- cosphi * y2 - sinphi * y3
x3 <- sinphi * y1
x4 <- cosphi * y1
R[, k, j] <- c(x1, x2, x3)/(1 - x4)
}
}
vs <- matrix(R, nrow = 3L, ncol = nu*nv)
tris <- cgalMeshes:::meshTopology(nu, nv, uperiodic, vperiodic)
tmesh3d(
vertices    = vs,
indices     = tris,
homogeneous = FALSE
)
}

If you run hMesh(gamma0, m) with m large enough, here is the mesh you will obtain (actually you have to close gamma0, that is to say you have to use rbind(gamma0, gamma0[1L, ])): Finally I did another animation. The Hopf torus whose profile curve is the equator of the unit sphere is nothing but an ordinary torus after the stereographic projection. Then, I scaled the $$x$$-coordinates of gamma0 continuously from a factor zero to a positive factor and I plotted the (stereographic projection of the) Hopf torus corresponding to each scaling. In this way the ordinary torus is smoothly transformed to the previous Hopf torus: The code:

h_ <- seq(0, 2, length.out = 60L) # scaling factors
open3d(windowRect = 50 + c(0, 0, 512, 512))
bg3d(rgb(54, 57, 64, maxColorValue = 255))
view3d(0, 0, zoom = 0.75)
for(i in seq_along(h_)) {
gamma <- gamma0
gamma[, 1L] <- h_[i] * gamma[, 1L]
# normalize so that the points are on the sphere
gamma <- gamma / sqrt(apply(gamma, 1L, crossprod))
gamma <- rbind(gamma, gamma[1L, ])
mesh <- hMesh(gamma, 500L)
mesh <- addNormals(mesh, angleWeighted = FALSE)
snapshot3d(sprintf("zzpic%03d.png", i), webshot = FALSE)
clear3d()
}
# mount the animation ####
library(gifski)
pngFiles <- Sys.glob("zzpic*.png")
gifski(
png_files = pngFiles,
gif_file = "HopfTorusCubicalConeToTorus.gif",
width  = 512,
height = 512,
delay = 1/15
)
file.remove(pngFiles)