Drawing a smooth implicit surface with misc3d

Question

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."?

Try adding `smooth=1000` to your `makeTriangles` statement, i.e. `drawScene.rgl(makeTriangles(surf, smooth=1000))` — G5W, Nov 13 '18 at 00:41
@G5W Oops sorry, I've just posted an answer with the option `smooth`, before seeing your comment...I will try `smooth=1000`. — Stéphane Laurent, Nov 13 '18 at 02:07
@G5W I've tried. It seems that `smooth=TRUE` is the same as `smooth=1000` when we use the `rgl` rendering. — Stéphane Laurent, Nov 13 '18 at 02:09
I tried another automatic approach: I used the suggestion from https://gamedev.stackexchange.com/a/145034 to approximate the normals at grid points by differencing and then interpolated for the triangle vertices (which always lie on a grid line). This seems to do a bit better than the `addNormals` approach based on the triangle normals, but still shows big artefacts. So I think if you want a nice smooth surface, you'll have to compute the true normals the way you did. — user2554330, Nov 16 '18 at 20:25

score 2 · Answer 1 · answered Nov 13 '18 at 01:41

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.

Stéphane Laurent · Answer 2 · 2022-06-03T06:13:36.527

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.

EDIT

The result is highly better with the rmarchingcubes package:

library(rmarchingcubes)
contour_shape <- contour3d(
  griddata = voxel, level = 0,
  x = x, y = y, z = z
)
library(rgl)
tmesh <- tmesh3d(
  vertices = t(contour_shape[["vertices"]]),
  indices = t(contour_shape[["triangles"]]),
  normals = contour_shape[["normals"]],
  homogeneous = FALSE
)
open3d(windowRect = c(50, 50, 562, 562))
view3d(zoom=0.8)
shade3d(tmesh, color = "darkred")

Yes, definitely preferable. – user2554330 Nov 13 '18 at 03:16 — user2554330, Nov 13 '18 at 03:16

Drawing a smooth implicit surface with misc3d

2 Answers2

EDIT