Calculating 3d contours

library(rmarchingcubes)

rmarchingcubes takes C++ code written by Thomas Lewiner for an efficient implementation of the marching cubes algorithm and ports into an R package. The key and only function exported is contour3d(), taking a 3-dimensional array of values and returning the calculated 3d mesh object fit to this data. A similar function is provided in the misc3d package with more flexibility for different inputs and outputs. The implementation here has two key advantages, firstly since the implementation is based on compiled C++ code the result should be considerably quicker, perhaps by orders of magnitude, secondly normals are additionally calculated and returned for each vertex making up the 3d contour.

A simple usage example

The following example goes through a simple example of generating some example input data and applying the algorithm. In this case values in the input data increase with the squared distance from the origin so the resulting 3d contour representing a given level is a sphere. In reality how you calculate your input data will vary depending on your intended usage.

The first job is to generate the required input data, a 3-dimensional array of values from which you wish to calculate a 3d contour mesh.

# Function to generate values decreasing in a sphere-like way
f <- function(coords) coords[1]^2 + coords[2]^2 + coords[3]^2

# Set grid coordinates at which to calculate values
x <- seq(-2,2,len = 20)
y <- seq(-2,2,len = 20)
z <- seq(-2,2,len = 20)

# Calculate values across grid coordinates
grid_coords <- expand.grid(x, y, z)
grid_values <- apply(grid_coords, 1, f)

# Convert to a 3d array
grid_array <- array(grid_values, dim = c(length(x), length(y), length(z)))

Now we can run contour3d() to apply the marching cubes algorithm and return mesh information representing a 3d contour at our desired level, in this case at a value of 4.

# Calculate 3d contour from the grid data at a contour level of value 4
contour_shape <- contour3d(
  griddata = grid_array, 
  level = 4,
  x = x,
  y = y,
  z = z
)

Finally once we have the mesh information we can view it using our desired method. In the example below I use the r3js package to view the 3d contour we have calculated.

# Optionally view the output using the r3js package
# devtools::install_github("shwilks/r3js")

# Setup plot object
data3js <- r3js::plot3js(
  x = x,
  y = y,
  z = z,
  type = "n"
)

# Add shape according to the calculated contours
data3js <- r3js::shape3js(
  data3js,
  vertices = contour_shape$vertices,
  faces = contour_shape$triangles,
  normals = contour_shape$normals,
  col = "red"
)

# View the plot
r3js::r3js(data3js)