How to clip an isosurface to a ball?

Question

Consider the Togliatti implicit surface. I want to clip it to the ball centered at the origin with radius 4.8. A solution, with the misc3d package, consists in using the mask argument of the computeContour3d function, which allows to use only the points satisfying x^2+y^2+z^2 < 4.8^2:

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))

But the borders of the resulting surface are irregular:

How to get regular, smooth borders?

Answer 1

Your solution is excellent for the problem you stated, because spherical coordinates are so natural for that boundary. However, here is a more general solution that would work for other smooth boundaries.

The idea is to allow input of a boundary function, and cull points when they are too large or too small. In your case it would be the squared distance from the origin, and you would want to cull points where the value is bigger than 4.8^2. But sometimes the triangles being drawn to make the smooth surface should only be partially culled: one point would be kept and two deleted, or two kept and one deleted. If you cull the whole triangle that leads to the jagged edges in your original plot.

To fix this, the points can be modified. If only one is supposed to be kept, then the other two points can be shrunk towards it until they lie on an approximation to the boundary. If two are supposed to be kept you want the shape to be a quadrilateral, so you would build that out of two triangles.

This function does that, assuming the input surf is the output of computeContour3d:

boundSurface <- function(surf, boundFn, bound = 0, greater = TRUE) {
  # Surf is n x 3:  each row is a point, triplets are triangles
  values <- matrix(boundFn(surf) - bound, 3)
  # values is (m = n/3) x 3:  each row is the boundFn value at one point
  # of a triangle
  if (!greater) 
    values <- -values
  keep <- values >= 0
  # counts is m vector counting number of points to keep in each triangle
  counts <- apply(keep, 2, sum)
  # result is initialized to an empty array
  result <- matrix(nrow = 0, ncol = 3)
  # singles is set to all the rows of surf where exactly one
  # point in the triangle is kept, say s x 3
  singles <- surf[rep(counts == 1, each = 3),]
  if (length(singles)) {
    # singleValues is a subset of values where only one vertex is kept
    singleValues <- values[, counts == 1]
    singleIndex <- 3*col(singleValues) + 1:3 - 3
    # good is the index of the vertex to keep, bad are those to fix
    good <- apply(singleValues, 2, function(col) which(col >= 0))
    bad <- apply(singleValues, 2, function(col) which(col < 0))
    for (j in 1:ncol(singleValues)) {
      goodval <- singleValues[good[j], j]
      for (i in 1:2) {
        badval <- singleValues[bad[i,j], j]
        alpha <- goodval/(goodval - badval)
        singles[singleIndex[bad[i,j], j], ] <- 
          (1-alpha)*singles[singleIndex[good[j], j],] +
             alpha *singles[singleIndex[bad[i,j], j],]
      }
    }
    result <- rbind(result, singles)
  }
  doubles <- surf[rep(counts == 2, each = 3),]
  if (length(doubles)) {
    # doubleValues is a subset of values where two vertices are kept
    doubleValues <- values[, counts == 2]
    doubleIndex <- 3*col(doubleValues) + 1:3 - 3
    doubles2 <- doubles
    # good is the index of the vertex to keep, bad are those to fix
    good <- apply(doubleValues, 2, function(col) which(col >= 0))
    bad <- apply(doubleValues, 2, function(col) which(col < 0))
    newvert <- matrix(NA, 2, 3)
    for (j in 1:ncol(doubleValues)) {
      badval <- doubleValues[bad[j], j]
      for (i in 1:2) {
        goodval <- doubleValues[good[i,j], j]
        alpha <- goodval/(goodval - badval)
        newvert[i,] <- 
          (1-alpha)*doubles[doubleIndex[good[i,j], j],] +
             alpha *doubles[doubleIndex[bad[j], j],]
      }
      doubles[doubleIndex[bad[j], j],] <- newvert[1,]
      doubles2[doubleIndex[good[1,j], j],] <- newvert[1,]
      doubles2[doubleIndex[bad[j], j],] <- newvert[2,]
    }
    result <- rbind(result, doubles, doubles2)
  }
  # Finally add all the rows of surf where the whole
  # triangle is kept
  rbind(result, surf[rep(counts == 3, each = 3),])
}

You would use it after computeContour3d and before makeTriangles, e.g.

fn <- function(x) { 
  apply(x^2, 1, sum)
}

drawScene.rgl(makeTriangles(boundSurface(surf, fn, bound = 4.8^2, 
                                         greater = FALSE), 
                            smooth = TRUE))

Here's the output I see:

It's not quite as good as yours, but it would work for many different boundary functions.

Edited to add: Version 0.100.26 of rgl now has a function clipMesh3d which incorporates these ideas.