R/get_value_from_ijk.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Value of the volume at a selection of DICOM indices
#' @description The \code{get.value.from.ijk} function calculates the value of 
#' a "volume" class object at DICOM indices i, j, k, whether they are 
#' integers or not.
#' @param ijk Vector or 3-column matrix of DICOM indices.
#' @param vol "volume" class object.
#' @param interpolate Boolean, default to \code{TRUE}. If \code{interpolate = TRUE}, a 
#' trilinear interpolation of the value of the voxels, relative to the values of 
#' adjacent voxels, is performed.
#' @param s.ijk Vector of 3 positive numbers greater than or equal to 1,
#' representing the new voxel ijk-sizes in which averaging is calculated.
#' @return Returns a vector of the values of the volume at the requested DICOM 
#' indices.
#' @seealso \link[espadon]{get.ijk.from.index}.
#' @examples
#' # loading of toy-patient objects (decrease dxyz and increase beam.nb for 
#' # better result)
#' step <- 4
#' patient <- toy.load.patient (modality = "rtdose", roi.name = "", 
#'                              dxyz = rep (step, 3), beam.nb = 3)
#' D <- patient$rtdose[[1]]
#' # isodose
#' Dmax <- max (D$vol3D.data, na.rm = TRUE)
#' Dmax
#' idx <- which (D$vol3D.data >= (Dmax -1) & D$vol3D.data <= (Dmax - 0.2))
#' ijk <- get.ijk.from.index (idx, D)
#' get.value.from.ijk (ijk, vol = D, interpolate = FALSE)

#' @export
#' @importFrom methods is
get.value.from.ijk <- function (ijk, vol, interpolate = TRUE, s.ijk = c(1,1,1))  {
  if (!is (vol, "volume")) stop ("vol should be a volume class object.")
  s_ijk <- abs(s.ijk)
  # s_ijk [s_ijk < 1] <- 1
  if (length(s_ijk) < 3) s_ijk <- c(s_ijk, rep(s_ijk[length(s_ijk) ],2))[1:3]
  # if (length(s.ijk)!=3 | !all(s_ijk[1:min(length(s.ijk),3)]==s.ijk[1:min(length(s.ijk),3)])) 
  #   warning("s.ijk is set to c(", paste(s_ijk, collapse = ", "),")")
  toC3M <- function (vect) {return(matrix(vect,ncol=3))}
  ijk <- toC3M (ijk)

  k.idx<- match(0:max(vol$k.idx),vol$k.idx)
  k.loc <- k.idx-1
  fna <- is.na(k.idx)
  k.idx[!fna] <- vol$k.idx
  k.loc[fna] <- max(vol$k.idx)+1
  k.idx[fna] <- max(vol$k.idx)+1
  vol3D <- as.numeric(vol$vol3D.data)
  # vol3D[is.na(vol3D)] <- NaN
  value <- .getvaluefromijkC (vol3D = vol3D,
                              interpolate = interpolate,
                              i = as.numeric(ijk[ ,1]),
                              j = as.numeric(ijk[ ,2]),
                              k = as.numeric(ijk[ ,3]),
                              k_idx = k.idx,
                              k_loc = k.loc, n_ijk=vol$n.ijk,
                              s_ijk = s_ijk)
  
  # value [is.nan(value)] <- NA
  return(value)
}