mesh3d from xyz points (with walls)

[mdsumner]

go crazy (needs expanding vertices for "vertices" colouring when walls present)

#' Create a mesh3d from xyz points
#'
#' Points are triangulated in x,y and then turned into mesh3d for the rgl package.
#' 
#' Walls cannot be made transparent - but plot3d(mesh, alpha = ) will work for global transparency). 
#' 
#' At the moment, colour will bleed down the walls because the vertices are shared ... WIP
#' @param xyz matrix of x,y,z values (points, no structure)
#' @param zwall optional walls added to the mesh, a numeric value for their 'base' (may be above or below z values)
#' @param ... ignored
#' @param meshColor see rgl::as.mesh3d ('vertices' for continuous colouring, 'faces' for discrete colours)
#' @param plot include a plot as a side-effect (TRUE by default)
#' @param walls the wall colour (defaults to grey, cannot be transparent
#'
#' @return mesh3d
#' @export
#'
#' @examples
#' points_to_mesh(cbind(quakes$long, quakes$lat, -quakes$depth), zwall = -500, meshColor = "faces")
points_to_mesh <- function(xyz, zwall = NULL,  ..., meshColor = "faces", plot = TRUE, walls = grey(0.5, alpha = 1)) {
  idx <- geometry::delaunayn(xyz[,1:2])
  mesh <- rgl::tmesh3d(t(cbind(xyz, h = 1)), t(idx), meshColor = meshColor)


  ## check it
  # anglr::mesh_plot(mesh)
  # par(xpd = NA)
  # lines(t(mesh$vb[1:2, t(cbind(boundary, NA))]), lwd = 5, col = "firebrick")
  wallcolor <- NULL
  if (meshColor == "faces") {
    cols <- colourvalues::color_values(colMeans(matrix(mesh$vb[3, mesh$it], 3L)))
    
  } else {
    # vertices
    cols <- colourvalues::color_values(mesh$vb[3, ])
  }
  if (!is.null(zwall)) {
    
    ## add walls
    edges <- t(cbind(mesh$it[1:2, ], mesh$it[2:3, ], mesh$it[c(3, 1), ]))
    native <- edges[,1] < edges[,2]
    ## by sorting each row
    edges1 <- cbind(pmin(edges[,1], edges[,2]), pmax(edges[,1], edges[,2]))
    ## and find duplicates, keep only those that occur once (they are mesh boundary)
    tib <- tibble::tibble(.vx0 = edges1[,1], .vx1 = edges1[,2])
    tib <- dplyr::ungroup(dplyr::mutate(dplyr::group_by(tib, .vx0, .vx1), n = dplyr::n() ))
    tib <- dplyr::filter(tib, n < 2)
    boundary <- as.matrix(tib[c(".vx0", ".vx1")])
    
    
    
    ## add the quads with a constant z
    nv <- rbind(mesh$vb[1:2, t(boundary[,2:1])], z = zwall, h = 1)
    ## or a relative z
    #nv <- rbind(mesh$vb[1:2, t(boundary[,2:1])], z = mesh$vb[3, rep(boundary[,1], each = 2)] - 30, h = 1)
    ncurr <- dim(mesh$vb)[2L]
    nnew <- nrow(boundary) * 2
    quads <- cbind(boundary, matrix(ncurr + seq_len(nnew), ncol = 2, byrow = TRUE))
    mesh$ib <- t(quads)            ## the quads are all new
    vb <- cbind(mesh$vb, nv)  ## but the geometry is part new
    if (meshColor == "faces") {
      nc <-  dim(boundary)[1L]
    } else {
      ## vertices
      nc <- dim(nv)[2L]
    }
    wallcolor <- rep(walls, nc)
    
    mesh$vb <- vb 
    
  }

 
  mesh$material$color <- c(cols, 
                         wallcolor) ##runif(ncol(mesh$ib))))
  if (plot) {
   rgl::open3d()
   rgl::plot3d(mesh)
  }
  mesh
}