Pfthb

The misc3d package provides a great implementation of the marching cubes algorithm, allowing to plot implicit surfaces.

For example, let's plot a Dupin cyclide:

a = 0.94; mu = 0.56; c = 0.34 # cyclide parameters
f <- function(x, y, z, a, c, mu) # implicit equation f(x,y,z)=0
 b <- sqrt(a^2-c^2)
 (x^2+y^2+z^2-mu^2+b^2)^2 - 4*(a*x-c*mu)^2 - 4*b^2*y^2

# define the "voxel"
nx <- 50; ny <- 50; nz <- 25
x <- seq(-c-mu-a, abs(mu-c)+a, length=nx) 
y <- seq(-mu-a, mu+a, length=ny)
z <- seq(-mu-c, mu+c, length=nz) 
g <- expand.grid(x=x, y=y, z=z)
voxel <- array(with(g, f(x,y,z,a,c,mu)), c(nx,ny,nz))
# plot the surface
library(misc3d)
surf <- computeContour3d(voxel, level=0, x=x, y=y, z=z)
drawScene.rgl(makeTriangles(surf))

Nice, except that the surface is not smooth.
The documentation of drawScene.rgl says: "Object-specific rendering features such as smoothing and material are controlled by setting in the objects." I don't know what does that mean. How to get a smooth surface?

I have a solution but not a straightforward one: this solution consists in building a mesh3d object from the output of computeContour3d, and to include the surface normals in this mesh3d.

The surface normals of an implicit surface defined by f(x,y,z)=0 are simply given by the gradient of f. It is not hard to derive the gradient for this example.

gradient <- function(xyz,a,c,mu)
 x <- xyz[1]; y <- xyz[2]; z <- xyz[3]
 b <- sqrt(a^2-c^2)
 c(
 2*(2*x)*(x^2+y^2+z^2-mu^2+b^2) - 8*a*(a*x-c*mu),
 2*(2*y)*(x^2+y^2+z^2-mu^2+b^2) - 8*b^2*y,
 2*(2*z)*(x^2+y^2+z^2-mu^2+b^2)
 )

Then the normals are computed as follows:

normals <- apply(surf, 1, function(xyz)
 gradient(xyz,a,c,mu)
)

Now we are ready to make the mesh3d object:

mesh <- list(vb = rbind(t(surf),1),
 it = matrix(1:nrow(surf), nrow=3),
 primitivetype = "triangle", 
 normals = rbind(-normals,1))
class(mesh) <- c("mesh3d", "shape3d")

And finally to plot it with rgl:

library(rgl)
shade3d(mesh, color="red")

Nice, the surface is smooth now.

But is there a more straightforward way to get a smooth surface, without building a mesh3d object? What do they mean in the documentation: "Object-specific rendering features such as smoothing and material are controlled by setting in the objects."?

I don't know what the documentation is suggesting. However, you can do it via a mesh object slightly more easily than you did (though the results aren't quite as nice), using the addNormals() function to calculate the normals automatically rather than by formula.

Here are the steps:

Compute the surface as you did.

Create the mesh without normals. This is basically what you did, but using tmesh3d():

mesh <- tmesh3d(t(surf), matrix(1:nrow(surf), nrow=3), homogeneous = FALSE)

Calculate which vertices are duplicates of which others:

verts <- apply(mesh$vb, 2, function(column) paste(column, collapse = " "))
firstcopy <- match(verts, verts)

Rewrite the indices to use the first copy. This is necessary, since the misc3d functions give a collection of disconnected triangles; we need to work out which are connected.

it <- as.numeric(mesh$it)
it <- firstcopy[it]
dim(it) <- dim(mesh$it)
mesh$it <- it

At this point, there are a lot of unused vertices in the mesh; if memory was a problem you might want to add a step to remove them. I'm going to skip that.

Add the normals

mesh <- addNormals(mesh)

Here are the before and after shots. Left is without normals, right is with them.

It's not quite as smooth as your solution using computed normals, but it's not always easy to find those.

There's an option smooth in the makeTriangles function:

drawScene.rgl(makeTriangles(surf, smooth=TRUE))

I think the result is equivalent to @user2554330's solution, but this is more straightforward.