R/fan_to_voxel.R

In espadon: Easy Study of Patient DICOM Data in Oncology

#' Indices of voxels crossed by a fan
#' @description The \code{fan.to.voxel} function computes the indices of voxels 
#' crossed by a fan. It is useful for retrieving voxel values and voxel indices 
#' of a volume (dose or ct) along the fan rays.
#' @param vol "volume" class object.
#' @param fan "fan" class object created by \link[espadon]{fan.sphere} for example.
# @param att Boolean. If \code{TRUE}, the cumulative sum  of the \code{vol} voxels 
# values along the ray path, times the distance crossed by the ray in these 
# voxels, is returned.
#' @param restrict Boolean. If \code{TRUE}, only the voxels with a value equal 
#' to \code{vol.value} are taken into account.
#' @param vol.value Value of the voxels taken into account, in case of \code{restrict = TRUE}
#' @return Returns a dataframe of 4 columns. Each line gives:
#' \itemize{
#' \item column "ray.index": the index (i.e. the row number) of the ray 
#' concerned in \code{fan$dxyz},
#' \item column "vol.index": the index of the voxel crossed in \code{vol$vol.3Ddata},
#' \item column "l.in": the distance between fan source (i.e. \code{fan$origin}) 
#' and the first face of the voxel crossed by the ray,
#' \item column "dl": the distance crossed by the ray in the voxel.
# \item if \code{att = TRUE} the cumulative sum of the voxels 
# values along the ray path, times dl (column "att").
#' }
#' If the rays do not cross any voxel, the dataframe has no row.
#' @seealso \link[espadon]{fan.beam}, \link[espadon]{fan.planar}, \link[espadon]{fan.sphere}. 
#' @export
#' @examples
#' vol <- vol.create (pt000 = c(1,10,10), dxyz = c (1 , 1, 1),
#'                    n.ijk = c(100, 100, 100)) 

#' fan.origin <- c (50,50,50)                         
#' fan <- fan.sphere (angle = 10, origin = fan.origin)

#' fan.voxel <- fan.to.voxel (vol = vol, fan = fan)
#' head (fan.voxel)
#' 
#' # Use of the 2nd column of fan.voxel to select voxels 
#' bin <- vol.copy (vol, modality = "binary")
#' bin$vol3D.data[] <- FALSE
#' bin$vol3D.data[fan.voxel[,2]] <- TRUE
#' bin$max.pixel <- TRUE
#' bin$min.pixel <- FALSE
#' display.kplane(bin, k=10)
 
fan.to.voxel <- function(vol, fan, restrict = FALSE, vol.value = 1) {

  if (!is(vol, "volume")) {
    stop("vol should be a volume class object.")
    # return (data.frame(ray.index=numeric(0), vol.index=numeric(0),
    #                    l.in=numeric(0), dl=numeric(0)))
  }
  
  if (!is(fan, "fan")) {
    stop("fan should be a fan class object.")
    # return (data.frame(ray.index=numeric(0), vol.index=numeric(0),
    #                    l.in=numeric(0), dl=numeric(0)))
  }
  
  if (fan$ref.pseudo != vol$ref.pseudo) {
    warning("vol ref.pseudo and fan ref.pseudo are different")
    # return (data.frame(ray.index=numeric(0), vol.index=numeric(0),
    #                    l.in=numeric(0), dl=numeric(0)))
  }
   
  #2D
  idx.c <- which(apply(abs(vol$xyz.from.ijk[1:3,1:3]),2,sum) == 0) 
  idx.r <-  which(apply(abs(vol$xyz.from.ijk[1:3,1:3]),1,sum) == 0)
  if (length(idx.c) > 0) {
    if (abs(vol$xyz0[1,idx.r]) > 1e-6) return(NULL)
    u <- vol$xyz.from.ijk 
    u[idx.r,idx.c] <- 1
    Mat <- solve(u)
    Mat[idx.r,idx.c] <- 0
  } else {#3D
    Mat <- solve(vol$xyz.from.ijk)}

  u_ijk <- (cbind(fan$xyz,0) %*% t(Mat))[,1:3]
  O_ijk <-  as.numeric((c(fan$origin,1) %*% t(Mat))[,1:3])
  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
  n_ijk  <- as.numeric(vol$n.ijk)
  
  p <- as.numeric(t(u_ijk))
  ncol <- 4
  # if (att) ncol <- 5
  df <- as.data.frame(matrix(.fantovoxelC(p, n_ijk, k_idx, 
                                          k_loc,O_ijk, 
                                          vol_data = as.numeric(vol$vol3D.data),
                                          att = FALSE,
                                          vol_value_flag = restrict,
                                          vol_value = as.numeric(vol.value)), 
                             ncol = ncol, byrow = TRUE))
  colnames(df) <- c("ray.index","vol.index","l.in","dl","att")[1:ncol]
  return(df)
}

.fan.to.voxel.with.att <- function(vol, fan, att,restrict = FALSE, vol.value = 1){

  if (!is(vol, "volume")) {
    stop("vol should be a volume class object.")
  }
  
  if (!is(fan, "fan")) {
    stop("fan should be a fan class object.")
  }
  
  if (fan$ref.pseudo != vol$ref.pseudo) {
    warning("vol ref.pseudo and fan ref.pseudo are different")
    # return (data.frame(ray.index=numeric(0), vol.index=numeric(0),
    #                    l.in=numeric(0), dl=numeric(0)))
  }
  
  Mat <- solve(vol$xyz.from.ijk) 
  u_ijk <- (cbind(fan$xyz,0) %*% t(Mat))[,1:3]
  O_ijk <-  as.numeric((c(fan$origin,1) %*% t(Mat))[,1:3])
  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
  n_ijk  <- as.numeric(vol$n.ijk)
  
  p <- as.numeric(t(u_ijk))
  ncol <- 4
  if (att) ncol <- 5
  df <- as.data.frame(matrix(.fantovoxelC(p, n_ijk, k_idx, 
                                          k_loc,O_ijk, 
                                          vol_data = as.numeric(vol$vol3D.data),
                                          att = att,
                                          vol_value_flag = restrict,
                                          vol_value = as.numeric(vol.value)), 
                             ncol = ncol, byrow = TRUE))
  colnames(df) <- c("ray.index","vol.index","l.in","dl","att")[1:ncol]
  return(df)
}