Slices of an implicit hypersurface with R
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
In a
previous post, I showed how to draw a 3D slice of a 4D hypersurface when a
parameterization of this hypersurface is available. Here we deal with
the case when an implicit equation of the hypersurface is available. For
the illustration, we again consider the
tiger. It is given by the implicit equation
(√x2+y2–R1)2+(√z2+w2–R2)2=r2.
# right-isoclinic rotation in 4D space # xi is the angle of 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] cbind( 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 ) }
So, in the implicit equation we fix w=w0 and we perform the rotation, taking arbitrary values for α and β:
f <- function(xyz, w0, xi, R1 = 3, R2 = 2, r = 1){ rxyzw <- rotate4d(pi/4, pi/4, xi, cbind(xyz, w0)) x <- rxyzw[, 1L] y <- rxyzw[, 2L] z <- rxyzw[, 3L] w <- rxyzw[, 4L] (sqrt(x^2+y^2) - R1)^2 + (sqrt(z^2+w^2) - R2)^2 - r^2 }
To plot the isosurface, we will use the rmarchingcubes package, not only for its speed, but also because it computes an excellent approximation of the per-vertex normals (it approximates the gradient of f). So taking a 150×150×150 grid is enough to get a smooth surface:
# make grid #### n <- 150L x <- seq(-5, 5, len = n) y <- seq(-5, 5, len = n) z <- seq(-5, 5, len = n) Grid <- expand.grid(X = x, Y = y, Z = z) # run the marching cubes #### library(rmarchingcubes) voxel <- array(f(Grid, w_0 = 0.3, xi = pi/3), dim = c(n, n, n)) cont <- contour3d(voxel, level = 0, x = x, y = y, z = z) # plot #### library(rgl) mesh <- tmesh3d( vertices = t(cont[["vertices"]]), indices = t(cont[["triangles"]]), normals = cont[["normals"]], homogeneous = FALSE ) open3d(windowRect = c(50, 50, 562, 562), zoom = 0.8) bg3d( sphere = FALSE, texture = "SpaceBackground.png", col = "white" ) shade3d(mesh, color = "maroon")
Now let’s make an animation by varying the angle of rotation ξ from 0 to π.
# vector of angles #### nframes <- 60L xi_ <- seq(0, pi, length.out = nframes) # open the 3D engine #### open3d( windowRect = c(50, 50, 562, 562), zoom = 0.85, userMatrix = rbind( c(0.93, -0.16, -0.33, 0), c(0.35, 0.66, 0.67, 0), c(0.11, -0.74, 0.67, 0), c( 0, 0, 0, 1) ) ) bg3d( sphere = FALSE, texture = "SpaceBackground.png", col = "white" ) # save the frames in png files #### for(i in 1L:nframes){ v <- array(f(Grid, w0 = 0.3, xi = xi_[i]), dim = c(n, n, n)) cont <- contour3d(v, level = 0, x = x, y = y, z = z) mesh <- tmesh3d( vertices = t(cont[["vertices"]]), indices = t(cont[["triangles"]]), normals = cont[["normals"]], homogeneous = FALSE ) shade3d(mesh, color = "maroon") snapshot3d(sprintf("zzpic%03d.png", i), webshot = FALSE) clear3d() } # make the animation with ImageMagick #### # option '-duplicate 1,-2-1' to get a forward-backward animation command <- paste0( "magick convert -dispose previous -delay 10 ", "-duplicate 1,-2-1 zzpic*.png tiger.gif" ) system(command)
A similar animation with a more complex surface can be found on my youtube channel.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.