R/vol_gradient.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Gradient of a volume
#' @description The \code{vol.gradient} function calculates the 3D gradient of a 
#' "volume" class object
#' @param vol "volume" class object.
#' @param alias Character string, \code{$alias} of the created object.
#' @param description Character string, describing the created object. If 
#' \code{description = NULL} (default value), it will be set to 
#' \code{paste (vol$object.alias, "gradient")}.
#' @return Returns a "volume" class object (see \link[espadon]{espadon.class} 
#' for class definitions), with the same grid and modality 
#' as \code{vol}, gradient of \code{vol}.
#' @examples
#' # loading of toy-patient objects (decrease dxyz and increase beam.nb for 
#' # better result)
#' step <- 3
#' pat <- toy.load.patient (modality = c ("ct"), dxyz = rep (step, 3))
#' CT.gradient <- vol.gradient (pat$ct[[1]])
#' 
#' display.plane (CT.gradient, view.type = "sagi", view.coord = 61, 
#'                interpolate = FALSE)

#' @export
vol.gradient <- function (vol,  alias = "", description = NULL) {
  
  if (is.null (vol)) return (NULL)
  if (!is (vol, "volume")) stop ("vol should be a volume class object.")
  if (is.null(vol$vol3D.data)) stop ("empty vol$vol3D.data.")
  
  if (is.null(description)) description <-  paste (vol$object.alias, "gradient")
  Gr <- vol.copy (vol, alias = alias, description = description)
  Gr$vol3D.data[] <- 0
  
  Grx <- Gry <- Grz <- Gr
  Gr$vol3D.data[] <- NA
  n <- vol$n.ijk[1]
  if (n>3) {
    Grx$vol3D.data[2:(n - 1),,] <- vol$vol3D.data[3:n,,] - vol$vol3D.data[1:(n - 2),,]
    Grx$vol3D.data <- Grx$vol3D.data / vol$dxyz[1]
    }
  n <- vol$n.ijk[2]
  if (n>3) {
    Gry$vol3D.data[,2:(n - 1),] <- vol$vol3D.data[,3:n,] - vol$vol3D.data[,1:(n - 2),]
    Gry$vol3D.data <- Gry$vol3D.data / vol$dxyz[2]
  }
  n <- vol$n.ijk[3]
  if (n>3) {
    Grz$vol3D.data[,,2:(n - 1)] <- vol$vol3D.data[,,3:n] - vol$vol3D.data[,,1:(n - 2)]
    Grz$vol3D.data <- Grz$vol3D.data / vol$dxyz[3]
  }
  Gr$vol3D.data <- sqrt (Grx$vol3D.data^2 + Gry$vol3D.data^2 + Grz$vol3D.data^2)
  # Gr$vol3D.data[c(1, Gr$n.ijk[1]),,] <- 0
  # Gr$vol3D.data[,c(1, Gr$n.ijk[2]),] <- 0
  # Gr$vol3D.data[,,c(1, Gr$n.ijk[3])] <- 0
  
  if (any(!is.na(Gr$vol3D.data))) {
    Gr$min.pixel <- min(Gr$vol3D.data, na.rm = TRUE)
    Gr$max.pixel <- max(Gr$vol3D.data, na.rm = TRUE)
  } else {
    Gr$min.pixel <- NA
    Gr$max.pixel <- NA
  }
  return (Gr)
  
}