Render mesh3d objects with rayshader and rayrender

Steps of the method from mesh3d to {rayshader}

Images below show the complete procedure on a small mesh with a low resolution:

• If needed, cut the 3d mesh with Morpho::cutMeshPlane()
• Project the mesh over a 2D raster with mesh2ray::mesh_to_raster()
• Transform the mesh as a matrix for {rayshader} with mesh2ray::stack_to_ray()
• Plot your matrix using {rayshader}

# Create simple mesh
simple_mesh_orig <- Rvcg::vcgIcosahedron()

# 1. Cut it with a plane with {Morpho}
v1 <- c(0, 0, -0.5)
v2 <- c(1, 0, -0.5)
v3 <- c(0, 1, -0.5)
simple_mesh <- Morpho::cutMeshPlane(simple_mesh_orig, v1, v2, v3)

# _plot3d
plot3d(simple_mesh, col = "cyan")
aspect3d(1, 1, 0.2)
plot3d(simple_mesh, type = "wire", add = TRUE)
plot3d(t(simple_mesh$vb), type = "s", radius = 0.1, add = TRUE)
view3d(theta = 20, phi = -40, zoom = 0.5, fov = 80)
# _snapshot
rgl.snapshot(filename = file.path(extraWD, "rayshader-simple-1-mesh.png"))
rgl::rgl.close()

# 2. Project on raster
simple_r <- mesh2ray::mesh_to_raster(simple_mesh, res = 26)
png(filename = file.path(extraWD, "rayshader-simple-2-raster.png"))
plot(simple_r)
dev.off()

# 3. raster to matrix
simple_ray <- mesh2ray::stack_to_ray(simple_r)

# 4. Rayshade matrix
zscale <- 0.4
ambmat <- ambient_shade(simple_ray$elevation, zscale = zscale)
raymat <- ray_shade(simple_ray$elevation, zscale = zscale, lambert = TRUE,
                    sunangle = 300)

ray_image <- simple_ray$elevation %>%
  sphere_shade(texture = "unicorn") %>%
  add_shadow(raymat, max_darken = 0.1) %>%
  add_shadow(ambmat, max_darken = 0.5)

# _plot 3d
ray_image %>% 
  plot_3d(simple_ray$elevation, zscale = zscale, windowsize = c(1000, 1000),
          soliddepth = -max(simple_ray$elevation, na.rm = TRUE)/zscale,
          theta = -30, phi = 50, zoom = 0.7, fov = 20
  )

# _snapshot
rgl.snapshot(filename = file.path(extraWD, "rayshader-simple-3-ray3d.png"))
rgl::rgl.close()

# Image montage
img_frames <- list.files(extraWD, pattern = "rayshader-simple", full.names = TRUE)
image_append(image_scale(image_read(img_frames), "400"))

Figure 1: Low resolution rayshading. a. mesh3d, b. rasterized mesh, c. rayshaded rasterized mesh

Transform a mesh3d as raster with {geometry}

The 3D human face can be projected in 2D. We can project it over a raster in the plane x/y. In the post with the 3D meristem, I used an inverse distance interpolation to estimate the height (z) of each cell of the projection raster. This calculation was quite heavy and I have hardly been able to use it on the human face mesh because it has a too high number of vertices (points). I tried different things with tools of my knowledge from {raster} package notably, but this was not enough. At least, I was able to produce something for my presentation “Everything but maps with spatial tools” at SatRday Paris on a low resolution version of the human face.
But I was not totally satisfied. I briefly discussed this on Twitter with Michael Sumner (mdsumner) who also like to play with 3D data and maps. A few hours later, he gave me a link to a script with a solution using geometry::tsearch. I enjoy so much how nice, friendly and welcoming is the R community !

Let’s build an empty raster of size 150*150, having the same x/y extent as our 3d face. The mesh3d object is a triangular mesh. In 2D, rach cell center of the empty raster falls in one of the triangles of the mesh. Using tsearch, you can find in which triangle of the mesh it falls. The function also returns the distance from the three angles of the triangle. With this information, you can interpolate linearly the z coordinates using a weighted mean (by the distance) of z of the three points of the triangle.

Note that this function is named mesh_to_raster() in my Github package {mesh2ray}

# If the orientation of the mesh please you,
# you can directly extract the xy vertices of the mesh triangles
# as a projection of the mesh on x/y
myMesh <- hum_mesh
triangle_xy <- t(myMesh$vb[1:2, ])
## Build a target raster with extent of the mesh
hum_mesh_r <- raster(extent(triangle_xy), nrows = 150, ncols = 150)
hum_mesh_r <- setValues(hum_mesh_r, NA_real_)
# Extract x/y coordinates of the raster
raster_xy <- sp::coordinates(hum_mesh_r)
## the magic barycentric index, the interpolation engine for the target grid vertices
## and where they fall in each triangle relative to triangle vertices
pid0 <- geometry::tsearch(x = triangle_xy[,1], y = triangle_xy[,2],
                          t = t(myMesh$it),
                          xi = raster_xy[,1], yi = raster_xy[, 2],
                          bary = TRUE)

# Determine raster cell center that really fall in a triangle of the mesh
ok <- !is.na(pid0$idx)
# Calculate linear interpolation of the z value using distances to vertices
hum_mesh_r[ok] <- colSums(matrix(myMesh$vb[3, myMesh$it[, pid0$idx[ok]]], nrow = 3) * t(pid0$p)[, ok])
# Plot the output raster
plot(hum_mesh_r)

Render the resulting raster using {rayshader}

That is easy to transform a raster to a matrix as used by {rayshader}, you simply need to transpose a as.matrix().

datamat <- t(as.matrix(hum_mesh_r))
# Rayshade raster
zscale <- 0.1
ambmat <- ambient_shade(datamat, zscale = zscale)
raymat <- ray_shade(datamat, zscale = zscale, lambert = TRUE,
                    sunangle = 45)

ray_image <- datamat %>%
  sphere_shade(texture = "unicorn") %>%
  add_shadow(raymat, max_darken = 0.1) %>%
  add_shadow(ambmat, max_darken = 0.5)

plot_map(ray_image)

If you want to overlay an image that has nothing to do with your raster, this may be a little tricky. In package {mesh2ray}, you can use function stack_to_ray() to resample an image to the dimension of the raster and get the matrix needed.
By the way, a {rayshader} output would not be a proper {rayshader} output without a gif…

r_img <- stack(file.path(extraWD, "ThinkR_logo_500px.png"))
hum_mesh_img <- stack_to_ray(hum_mesh_r, r_img)

# Rayshade raster
zscale <- 0.25
ambmat <- ambient_shade(hum_mesh_img$elevation, zscale = zscale)
raymat <- ray_shade(hum_mesh_img$elevation, zscale = zscale, lambert = TRUE)

ray_image_overlay <- hum_mesh_img$elevation %>%
  sphere_shade(texture = "imhof4") %>%
  add_overlay(hum_mesh_img$overlay, alphalayer = 0.99) %>% # Note overlay
  add_shadow(raymat, max_darken = 0.01) %>%
  add_shadow(ambmat, max_darken = 0.5) # %>% 
# add_water(watermap, color = "#57B6FF") 

# transform as gif
path_gif <- file.path(extraWD, "face.gif")
if (!file.exists(path_gif)) {
  # Add some movement
  th <- seq(-50, 50, length.out = 20)
  ph <- seq(80, 90, length.out = 20)
  move <- data.frame(theta = c(th, rev(th)), phi = c(ph, rev(ph)))
  
  n_frames <- nrow(move)
  dir.create(file.path(extraWD, "face_gif"), showWarnings = FALSE)
  img_frames <- file.path(
    extraWD, "face_gif",
    paste0("face", formatC(seq_len(n_frames), width = 3, flag = "0"), ".png"))
  
  # create all frames
  for (i in seq_len(n_frames)) {
    ray_image_overlay %>% 
      plot_3d(hum_mesh_img$elevation, zscale = zscale, windowsize = c(1000, 1000),
              soliddepth = -max(hum_mesh_img$elevation, na.rm = TRUE)/zscale,
              theta = move$theta[i], phi = move$phi[i], zoom = 0.30, fov = 20
      )
    
    rgl::snapshot3d(img_frames[i])
    rgl::clear3d()
  }
  rgl::rgl.close()
  
  magick::image_write_gif(magick::image_read(img_frames), 
                          path = path_gif, 
                          delay = 6/n_frames)
}
file.copy(file.path(extraWD, "face.gif"), 
          file.path(glue("../../static{StaticImgWD}/face.gif")), 
          overwrite = TRUE)

Steps of the method from mesh3d to {rayrender}

Images below show the complete procedure on a small mesh with a low resolution:

• Uniformely re-sampled the 3d mesh with Rvcg::vcgUniformRemesh
• Extract information of the uniform mesh with mesh2ray::mesh_to_cubes (that includes vcgUniformRemesh)
• Add cubes to the scene with mesh2ray::add_cubes_to_scene
• Plot the scene using {rayrender}

# Create simple mesh
simple_mesh <- Rvcg::vcgIcosahedron()

# _plot3d
plot3d(simple_mesh, col = "cyan")
aspect3d("iso")
plot3d(simple_mesh, type = "wire", add = TRUE)
plot3d(t(simple_mesh$vb), type = "s", radius = 0.1, add = TRUE)
view3d(theta = 0, phi = 10, zoom = 0.6, fov = 80)
# _snapshot
rgl.snapshot(filename = file.path(extraWD, "rayrender-simple-1-mesh.png"))
rgl::rgl.close()

# 1. Uniformely resample mesh
# _vcgUniformRemesh is included in mesh_to_cubes()
simple_resample <- vcgUniformRemesh(simple_mesh, voxelSize = 0.5, multiSample = TRUE,
                                discretize = TRUE)

plot3d(simple_resample, col = "cyan")
aspect3d("iso")
plot3d(simple_resample, type = "wire", add = TRUE)
plot3d(t(simple_resample$vb), type = "s", radius = 0.1, add = TRUE)
view3d(theta = 0, phi = 10, zoom = 0.6, fov = 80)
# _snapshot
rgl.snapshot(filename = file.path(extraWD, "rayrender-simple-2-regularmesh.png"))
rgl::rgl.close()

# 2. Calculate cube positions in a specific scene
# _vcgUniformRemesh is included in mesh_to_cubes()
simple_cubes <- mesh2ray::mesh_to_cubes(simple_mesh, voxelSize = 0.5, scene_dim = c(80, 475))

# 3. Draw with rayrender
# _scene
scene <- generate_cornell(lightintensity = 10)
# _add cubes on scene
scene <- mesh2ray::add_cubes_to_scene(scene, cubes = simple_cubes, 
                                      material = dielectric(color = "green"))
# _draw scene /!\ long calculation /!\
if (!file.exists(file.path(extraWD, "rayrender-simple-3-scene.png"))) {
  options(cores = 4)
  render_scene(scene, lookfrom = c(278, 278, -800) ,
               lookat = c(278, 278, 0), fov = 40, ambient_light = FALSE,
               samples = 500, parallel = TRUE, clamp_value = 5,
               filename =  file.path(extraWD, "rayrender-simple-3-scene.png"))
}

# Image montage
img_frames <- list.files(extraWD, pattern = "rayrender-simple", full.names = TRUE)
image_append(image_scale(image_read(img_frames), "400"))

Figure 2: Low resolution rayrendering. a. mesh3d, b. regularly sampled mesh, c. rayrendered regularly sampled mesh

Apply voxelisation on the human face and defines colours

The steps are the same than the small example above. We are using a higher resolution to get a final form much closer to the original.
However, I wanted to add some colours to some of the cubes and to do so, I did some manual drawing on 2D plots. After a few spatial points and polygons manipulations using {sf}, I can have a colored version of the face in 2D.

# voxelisation
data(humface)
human_cubes <- mesh2ray::mesh_to_cubes(humface, voxelSize = 0.5, scene_dim = c(20, 535))

# Define colours areas manually
polygon_contours_file <- file.path(extraWD, "polygon_contours.rds")
if (!file.exists(polygon_contours_file)) {
  # Manually draw polygons
  plot(human_cubes$cube_center$x, human_cubes$cube_center$y, pch = 20, cex = 0.5)

  eye1 <- raster::drawPoly(col = "blue")
  eye2 <- raster::drawPoly(col = "blue")
  mouth <- raster::drawPoly(col = "red")
  front <- raster::drawPoly(col = "green")
  
  polygon_contours <- list(eye1 = eye1, eye2 = eye2, mouth = mouth, front = front)
  # Transform as sf polygon
  polygon_contours_sf <- polygon_contours %>% 
  purrr::map2(names(polygon_contours), 
              ~sf::st_as_sf(.x) %>% mutate(name = .y)) %>% 
  do.call("rbind", .)
  
  readr::write_rds(polygon_contours_sf, polygon_contours_file)
} else {
  polygon_contours_sf <- readr::read_rds(polygon_contours_file)
}

# Transform cubes as sf for intersection with contours
human_colours <- human_cubes$cube_center %>% 
  sf::st_as_sf(coords = c("x", "y"), remove = FALSE) %>%
  sf::st_join(polygon_contours_sf) %>% 
  mutate(colours = case_when(
    name %in% c("eye1", "eye2") ~ "blue",
    name == "mouth" ~ "red",
    name == "front" ~ "ivory",
    TRUE ~ "grey"
  )) %>% 
  as_tibble()

# Explore the 3D points and color
plot3d(t(human_cubes$mesh_resample$vb), type = "s", radius = 1,
       col = human_colours$colours)
aspect3d("iso")
view3d(theta = 0, phi = 10, zoom = 0.6, fov = 80)
# _snapshot
rgl.snapshot(filename = file.path(tmpImgWD, "figure-html/humcolours.png"))
rgl::rgl.close()

Define material for each cube

I can create a tibble containing a material for each of the cubes to draw.

# All types of cubes
colours_cubes <- bind_rows(
  lambertian(color = "ivory"),
  dielectric(color = "#2E7BE6"), 
  dielectric(color = "#2E7BE6"), 
  dielectric(color = "#DB3B3B"),
  dielectric()
) %>% 
  mutate(name = c(NA, "eye1", "eye2", "mouth", "front"))

# Join with table of colours
human_materials <- dplyr::select(human_colours, name) %>% 
  left_join(colours_cubes, by = "name") %>% 
  dplyr::select(-name)

human_materials

## # A tibble: 245,975 x 12
##    type  properties checkercolor noise noisephase noiseintensity noisecolor
##    <chr> <list>     <list>       <dbl>      <dbl>          <dbl> <list>    
##  1 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  2 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  3 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  4 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  5 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  6 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  7 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  8 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
##  9 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
## 10 lamb… <dbl [3]>  <dbl [2]>        0          0             10 <dbl [3]> 
## # … with 245,965 more rows, and 5 more variables: image <list>,
## #   lightintensity <lgl>, fog <lgl>, fogdensity <dbl>, implicit_sample <lgl>

Draw the complete scene

The voxelised human_cubes object needs to be added to the empty scene using mesh2ray::add_cubes_to_scene with the tibble of materials created. Drawing the scene requires some time… Be patient… You can parallelise the calculations to reduce time.

# 3. Draw with rayrender
if (!file.exists(file.path(extraWD, "rayrender-human-scene.png"))) {
  # _scene
  scene <- generate_cornell(lightintensity = 20)
  # _add cubes on scene with material data.frame /!\ A little long /!\
  scene <- mesh2ray::add_cubes_to_scene(scene, cubes = human_cubes, 
                                        material = human_materials)
  
  # _draw scene /!\ long calculation /!\
  options(cores = 4)
  render_scene(scene, width = 600, height = 600,
               lookfrom = c(278, 278, -800) ,
               lookat = c(278, 278, 0), fov = 40, ambient_light = FALSE,
               samples = 500, parallel = TRUE, clamp_value = 5,
               filename =  file.path(tmpImgWD, "figure-html/rayrender-human-scene.png"))
}

Tada ! Now, it’s your turn. Although this two {rayshader} and {rayrender} packages are young, there are already many possibilities. I am sure there will be a lot of great inventions around them. And always remember to take some time to have fun ! R is also here for that !