terrainmesh to mesh3d with normals

as.mesh3d.terrainmeshr.R

library(terrainmeshr)
library(rgl)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Create triangles from heigh-map
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tris <- terrainmeshr::triangulate_matrix(
  volcano, maxError = 30, verbose = TRUE, y_up = FALSE
)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Unique vertex numbering
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
hash <- paste(tris[,1], tris[,2], tris[,3], sep = "-")
vidx <- as.integer(as.factor(hash))

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# calculate:
#    vb - the unique list of vertices (in order) 
#    it - the vertex indices which make up each face
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tris <- cbind(tris, vidx)

idx <- !duplicated(tris[,5])
vb  <- tris[idx,]
vb  <- vb[order(vb[,5]),]
vb  <- t(vb[,1:3])
vb  <- rbind(vb, 1)

it <- matrix(tris[,5], nrow = 3)


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Create a mesh3d object
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
object <- list(
  vb        = vb, 
  it        = it, 
  material  = list(), 
  normals   = NULL, 
  texcoords = NULL,
  meshColor = 'vertices'
)
class(object) <- c("mesh3d", "shape3d")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extract vertex normals as a matrix of line segments
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
get_normals_as_lines <- function(object, len = 10) {
  start <- object$vb[-4,]
  delta <- object$normals[-4,]
  end   <- start + delta * len
  lines <- matrix(rbind(start, end), nrow = 3)
  lines <- t(lines)
  lines
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Calculate normals with `rgl`
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
object <- rgl::addNormals(object)
normals <- get_normals_as_lines(object)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Shade and light it
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
open3d(windowRect = c(0,0, 800, 800))
shade3d(object, lit=TRUE, color='grey90', smooth=TRUE)
aspect3d(1, 1, 0.6)
rgl.viewpoint(-10, -50)

segments3d(normals, color = 'red', lwd = 3)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Calculate normals with Rvcg
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
object <- Rvcg:::vcgUpdateNormals(object, type = 0)
normals <- get_normals_as_lines(object)


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Shade and light it
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
open3d(windowRect = c(0, 0, 800, 800))
shade3d(object, lit=TRUE, color='grey90', smooth=TRUE)
aspect3d(1, 1, 0.6)
rgl.viewpoint(-10, -50)

segments3d(normals, color = 'red', lwd = 3)