R/nesting_bin.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Restrict volume to a binary selection
#' @description The \code{nesting.bin} function restricts a "volume" class 
#' object to the rectangular parallelepiped circumscribed to the selected voxels.
#' @param obj "volume" class object, containing data to restrict.
#' @param sel.bin "volume" class object, of "binary" modality, specifying the selected voxels.
#' @param xyz.margin Numeric vector of length 3, by default set to \code{c (0, 0, 0)}.
#' See details.
#' @param obj.restrict Boolean. Used if \code{obj} is of class "volume". If 
#' \code{obj.restrict = TRUE}, the rectangular parallelepiped circumscribed to 
#' the selected voxels, enlarged by xyz.margin cannot exceed the initial volume.
#' @param alias Character string, \code{$alias} of the created object.
#' @param description Character string, describing the the created object. 
#' If \code{description = NULL}, it will be 
#' \code{paste (obj$description,"restricted to", sel.bin$description)}.
#' @param ... Argument such as T.MAT, or for deprecated arguments
#' @return Returns a "volume" class object, in which 3D volume is limited to the
#' rectangular parallelepiped circumscribed to the voxels selected by \code{sel.bin}, increased by the
#' requested margins.

#' @details If \code{obj} is of class "volume", \code{xyz.margin} represents the 
#' distances in mm to be added to the x, y and z directions of the rectangular 
#' parallelepiped circumscribed to the voxels selected in \code{sel.bin}, in the 
#' \code{obj} frame of reference.
#' 
#' @details If \code{obj} is of class “mesh”, \code{sel.bin} will undergo a 
#' margin expansion \code{xyz.margin} before the mesh points are selected..

#' @seealso \link[espadon]{add.margin}, \link[espadon]{nesting.cube} and  
#' \link[espadon]{nesting.roi}.
#' @examples
#' # loading of toy-patient objects (decrease dxyz for  better result)
#' step <- 4
#' patient <- toy.load.patient(modality = c("ct", "rtstruct"), 
#'                             roi.name = "brain", dxyz = rep(step, 3))
#' CT <- patient$ct[[1]]
#' b <- bin.from.vol(CT, min = 0, max =200)
#' 
#' CT.restricted <- nesting.bin(CT, b, xyz.margin = rep (step, 3))
#' display.plane(bottom = CT.restricted, top = b, view.type = "sagi",
#'               bottom.col = pal.RVV(1000),
#'               bottom.breaks = seq(-1000, 1000, length.out = 1001),
#'               bg = "#00ff00",  interpolate  = FALSE)
#' @export              
nesting.bin <- function(obj, sel.bin, alias = "", description = NULL, 
                        xyz.margin = c(0, 0, 0), obj.restrict = FALSE, ...) {
  
  passed <- names(as.list(match.call())[-1])
  args <- list(...)
  if (!("obj" %in% passed)) {
    if (is.null(args[['vol']])) stop('argument "obj" is missing, with no default')
    obj <- args[['vol']]
  }

  if (!is.null(args[['vol.restrict']])) obj.restrict <- args[['vol.restrict']]
  eps <- args[["eps"]]
  if (is.null(eps)) eps <- 1e-9
  
  T.MAT <- args[['T.MAT']]
  
  if (!is(sel.bin, "volume")) stop("sel.bin should be a volume class object.")
 
  if (is.null(sel.bin$vol3D.data)) stop("empty sel.bin$vol3D.data.")
  
  if (is.null(description)) description <- paste(obj$description, "restricted to", 
                                                sel.bin$description)
  
  alias.sel <- sel.bin$object.alias
  sel.bin <- vol.in.new.ref(sel.bin,obj$ref.pseudo, T.MAT)
  sel.bin$object.alias <- alias.sel
  
  if (is(obj, "volume")) {
  #-----------------------
    if (sel.bin$modality != "binary" & sel.bin$modality != "weight") 
      stop("sel.bin must be of modality binary or weight.")
    pt <- get.xyz.from.index(which(sel.bin$vol3D.data > 0), sel.bin)
    if (is.null(pt)) stop("sel.bin has no selection zone")
    ext.pt <- data.frame(min = apply(pt,2,min), max = apply(pt,2,max))
    # ext.pt <- get.extreme.pt(sel.bin, pt.min = eps, min=0.5)
    ext.pt[1,] <- ext.pt[1,] + abs(xyz.margin[1])*c(-1,1)
    ext.pt[2,] <- ext.pt[2,] + abs(xyz.margin[2])*c(-1,1)
    ext.pt[3,] <- ext.pt[3,] + abs(xyz.margin[3])*c(-1,1)
    
    ext.pt.obj <- get.extreme.pt(obj)
    if (obj.restrict) {
      ext.pt[1,] <- .range.intersection(ext.pt[1,],ext.pt.obj[1,], egal = c(TRUE,TRUE))
      ext.pt[2,] <- .range.intersection(ext.pt[2,],ext.pt.obj[2,], egal = c(TRUE,TRUE))
      ext.pt[3,] <- .range.intersection(ext.pt[3,],ext.pt.obj[3,], egal = c(TRUE,TRUE))
    }
    if (any(is.na(ext.pt))) stop("obj and the selection of sel.bin have no common zones")
    obj_ <- nesting.cube(obj,pt.min = ext.pt[,1], pt.max = ext.pt[,2],
                         alias = alias, description = description)
    
  } else if (is(obj, "mesh")) {
    if (sel.bin$modality != "binary") stop("sel.bin must be of binary modality.")
    #on aggrandit sel.bin de la marge: 
  
    if (any(abs(xyz.margin) != 0)) sel.bin <- bin.dilation(sel.bin, radius = abs(xyz.margin))
    sel.bin$object.alias <- alias.sel
    
    obj_ <- obj
    obj_$description <- description
    obj_$object.alias <- alias
    if (obj_$nb.faces == 0)  return(obj_)
    
    if ((obj$ref.pseudo != sel.bin$ref.pseudo) & is.null(T.MAT)) 
      stop("obj and sel.bin have different ref.pseudo. Set T.MAT argument")
    f <- get.value.from.xyz(t(obj_$mesh$vb[1:3,]), sel.bin, T.MAT = T.MAT) >= 0.5
    f[is.na(f)] <- FALSE
    to.suppress <- which(!f)
    obj_$mesh$normals <- obj_$mesh$normals[,f]
    obj_$mesh$vb <- obj_$mesh$vb[,f]
    to.keep <- which(f)
    face.keep <- is.na(match(obj_$mesh$it[1,],to.suppress)) & 
      is.na(match(obj_$mesh$it[2,],to.suppress)) & 
      is.na(match(obj_$mesh$it[3,],to.suppress))
    obj_$mesh$it <- obj_$mesh$it[ ,face.keep]
    obj_$mesh$it[1,] <- match(obj_$mesh$it[1,],to.keep)
    obj_$mesh$it[2,] <- match(obj_$mesh$it[2,],to.keep)
    obj_$mesh$it[3,] <- match(obj_$mesh$it[3,],to.keep)
    if (!is.null(obj_$mesh$remface)) obj_$mesh$remface <- obj$mesh$remface[face.keep]
   
  }
  if (alias == "") return(obj_)
  return(.set.ref.obj(obj_,list(obj,sel.bin), add = FALSE))
}