R/mesh_from_bin.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Creation of a mesh according to a binary volume
#' @description The \code{mesh.from.bin} function creates a \code{mesh} class 
#' object from a volume object of "binary" modality.
#' @param bin "volume" class object of "binary" modality.
#' @param alias Character string, \code{$alias} of the mesh defining 
#' the \code{$alias} of the created object.
#' @param tol Tolerance in mm, applied for mesh simplification. See \link[Rvcg]{vcgClean}. 
#' The default value, equal to half the smallest voxel edge, limits meshing errors.
#' @param smooth.iteration Number of iterations applied in the smoothing algorithm. 
#' See \link[Rvcg]{vcgSmooth}. 
#' @param smooth.type character: select smoothing algorithm. Available are "taubin", 
#' "laplace", "HClaplace", "fujiLaplace", "angWeight" (and any sensible abbreviations). 
#' By default, set to "taubin". See \link[Rvcg]{vcgSmooth}.
#' @param smooth.lambda numeric: parameter for Taubin smooth. See \link[Rvcg]{vcgSmooth}.
#' @param smooth.mu numeric: parameter for Taubin smooth. See \link[Rvcg]{vcgSmooth}.
#' @param smooth.delta numeric: parameter for Scale dependent laplacian smoothing 
#' (see reference below).and maximum allowed angle (in radians) for deviation 
#' between normals Laplacian (surface preserving). See \link[Rvcg]{vcgSmooth}.
#' @param verbose Boolean, by default set to \code{FALSE}. Allows you to inhibit comments.
#' @seealso \link[Rvcg]{vcgSmooth} 
#' @return Returns a "mesh" class object. This is a list including the
#' following 6 elements:
#' \itemize{
#' \item \code{$patient}: set to \code{bin$patient}.
#' \item \code{$patient.bd}: set to \code{bin$patient.bd}.
#' \item \code{$patient.name}: set to \code{bin$patient.name}.
#' \item \code{$patient.sex}: set to \code{bin$patient.sex}.
#' \item \code{$file.basename}: set to "".
#' \item \code{$file.dirname}: set to "".
#' \item \code{$object.name}: set to "".
#' \item \code{$object.alias}: set to the \code{alias} argument of the function.
#' \item \code{$frame.of.reference}: set to \code{bin$frame.of.reference}.
#' \item \code{$ref.pseudo} : set to \code{bin$ref.pseudo}.
#' \item \code{$modality} : set to \code{"mesh"}.
#' \item \code{$description}: By default, set to \code{paste (bin$object.alias, "mesh")}.
#' \item \code{$creation.date}: set to \code{Sys.Date}.
#' \item \code{$nb.faces}: set to the number of faces of the mesh.
#' \item \code{$mesh}: list of 3 elements defining the mesh:
#' \tabular{rl}{
#' \tab - \code{$vb}: array made up of the generalized coordinates (x, y, z, 1) of the vertices of the triangles.\cr 
#'  \tab There are as many columns as there are vertices.\cr
#' \tab - \code{$it}: array of the 3 indices of the vertices forming a triangle, arranged by column.\cr
#'  \tab There are as many columns as there are triangles in the mesh.\cr
#' \tab - \code{$normals}: array made up of the generalized coordinates (x, y, z, 1) of the normal vectors of each triangle.\cr
#'  \tab There are as many columns as there are vertices.\cr
#' }
#' }
#' @note To compute the mesh, all NA voxels of the binary volume \code{bin} are 
#' set to FALSE. If all voxels are equal to FALSE, the function returns the code \code{NULL}.
#' @examples
#' # loading of toy-patient objects (decrease dxyz for better result)
#' step <- 4
#' patient <- toy.load.patient(modality = c("ct", "rtstruct"), roi.name = "",
#'                             dxyz = rep(step, 3))
#' CT <- patient$ct[[1]]
#' S <- patient$rtstruct[[1]]
#' 
#' # creation of the patient mesh
#' bin <- bin.from.roi(CT, struct = S, roi.name = "patient", verbose = FALSE)
#' mesh.patient <- mesh.from.bin(bin, alias = "patient", verbose = FALSE)
#' str(mesh.patient)

#' @importFrom Rvcg vcgClean vcgIsolated vcgSmooth
#' @importFrom misc3d contour3d
#' @importFrom rgl tmesh3d
#' @importFrom methods is
#' @export
mesh.from.bin <- function(bin, alias="", tol = min(abs(bin$dxyz))/2, 
                          smooth.iteration = 10, 
                          smooth.type = c("taubin", "laplace", "HClaplace", 
                                          "fujiLaplace", "angWeight",
                                          "surfPreserveLaplace"),
                          smooth.lambda = 0.5,
                          smooth.mu = -0.53,
                          smooth.delta = 0.1,
                          verbose = FALSE) {
  
  if (!is(bin, "volume")) stop("bin should be a volume class object.")
  if ((bin$modality != "binary") & (bin$modality != "weight")) stop("bin must be modality binary or weight.")

  vol3D <- bin$vol3D.data
  vol3D[,,1] <- FALSE
  vol3D[,,dim(vol3D)[3]] <- FALSE
  vol3D[is.na(vol3D)] <- FALSE
  
  if (all(vol3D == FALSE)) return(NULL)
  mesh_ <- list()
  
  mesh_$patient <- bin$patient
  mesh_$patient.name <- bin$patient.name
  mesh_$patient.bd <- bin$patient.bd
  mesh_$patient.sex <- bin$patient.sex
  
  mesh_$file.basename <- ""
  mesh_$file.dirname <- ""
  mesh_$object.name <- alias
  mesh_$object.alias <- alias
  mesh_$frame.of.reference <- bin$frame.of.reference
  mesh_$ref.pseudo <- bin$ref.pseudo
  mesh_$modality <- "mesh"
  mesh_$description <- paste(bin$object.alias, "mesh")
  mesh_$creation.date <- format(Sys.Date(), "%Y%m%d")

  mesh_$nb.faces <- 0
  
  mesh <- contour3d(vol3D,
                    x = 1:bin$n.ijk[1] - 1, y = 1:bin$n.ijk[2] - 1, z = bin$k.idx,
                    .5, fill = TRUE, smooth = TRUE, draw = FALSE)
  
  
  vertices <- matrix(as.vector(t(cbind(mesh$v1, mesh$v2, mesh$v3))), ncol = 3, byrow = TRUE)
  
  vertices <- cbind(vertices, rep(1, nrow(vertices)))
  N.triangles <- nrow(vertices) / 3
  indices <- cbind(3 * (1:N.triangles) - 2, 3 * (1:N.triangles) - 1, 3 * (1:N.triangles))
  mesh <- tmesh3d(t(vertices), t(indices))
  
  if (verbose) cat("mesh.from.bin INFO : Basic mesh cleanings\n")
  mesh <- vcgIsolated(mesh, silent = !verbose)
  mesh <- vcgClean(mesh, sel = c(0:4,7), tol = tol, iterate = TRUE, silent = !verbose)

  if (smooth.iteration > 0) {
    if (verbose) cat("mesh.from.bin INFO : mesh smoothing\n")
    smooth.type <- smooth.type[1]
    mesh <- vcgSmooth(mesh, iteration = smooth.iteration, type = smooth.type,
                      lambda = smooth.lambda, mu = smooth.mu,
                      delta = smooth.delta)
  }

  mesh$vb <- bin$xyz.from.ijk %*% mesh$vb
  mesh$remvert <- NULL
  mesh_$mesh <- mesh
  mesh_$nb.faces <- dim(mesh$it)[2]
  class(mesh_) <- "mesh"
  if (alias == "") return(mesh_)
  return(.set.ref.obj(mesh_,list(bin)))
  
}