library(misc3d)

# Togliatti surface equation: f(x,y,z) = 0
f <- function(x,y,z){
  w <- 1
  64*(x-w)*
    (x^4-4*x^3*w-10*x^2*y^2-4*x^2*w^2+16*x*w^3-20*x*y^2*w+5*y^4+16*w^4-20*y^2*w^2) -
    5*sqrt(5-sqrt(5))*(2*z-sqrt(5-sqrt(5))*w)*(4*(x^2+y^2-z^2)+(1+3*sqrt(5))*w^2)^2
}

# make grid
nx <- 220; ny <- 220; nz <- 220
x <- seq(-5, 5, length=nx)
y <- seq(-5, 5, length=ny)
z <- seq(-4, 4, length=nz)
g <- expand.grid(x=x, y=y, z=z)

# calculate voxel
voxel <- array(with(g, f(x,y,z)), dim = c(nx,ny,nz))

# mask: keep points satisfying x^2+y^2+z^2 < 4.8^2, in order to
#       clip the surface to the ball of radius 4.8
mask <- array(with(g, x^2+y^2+z^2 < 4.8^2), dim = c(nx,ny,nz))

# compute isosurface
surf <- computeContour3d(voxel, maxvol=max(voxel), level=0, mask=mask, x=x, y=y, z=z)

# draw isosurface
drawScene.rgl(makeTriangles(surf, smooth=TRUE))