R/amoment.R

In neuroconductor-devel/ichseg: Intracerebral Hemorrhage Segmentation of X-Ray Computed Tomography (CT) Images

#' @title Calculate Moments of Neighborhood
#'
#' @description This function calculates Local Moments (mean, standard deviation, skew) for an array.
#' @param image input image
#' @param window window (in width) for the neighborhood
#' @param nvoxels window (in voxels) for the neighborhood 1 results in a 3x3 cube
#' @param moment vector moments taken (1- mean, 2-sd, 3-skew)
#' @param mask array or object of class nifti of same size as image
#' @param only.mask Should objects outside the mask (i.e. zeros) be
#' counted the moment?  Default is FALSE so edges are weighted to 0
#' @param center vector of indicator of central moment.
#' if TRUE mean image is subtracted.  Same length as moment
#' @param invert Standardize the values by the power: 1/moment
#' @param mean_image mean image to be subtracted. If not supplied, and central = TRUE, local_moment_edge is run with mom = 1
#' @param na.rm remove NAs from the moment image calculation
#' @param remask set areas outside of mask to 0
#' @param ... Arguments passed to \code{\link{get_neighbors}}
#' @importFrom magic ashift
#' @importFrom neurobase niftiarr datatyper
#' @export
#' @return List of arrays same lenght as moment
#' @examples
#' x = array(rnorm(1000), dim=c(10, 10, 10))
#' mask = abs(x) < 1
#' mean.x = local_moment(x, nvoxels=1, moment = 1, mask=mask,
#' center = FALSE,
#' remask = FALSE)[[1]]
#' var.x = local_moment(x, nvoxels=1, moment = 2, mask=mask, center = TRUE,
#' mean_image = mean.x, remask=FALSE)[[1]]
#'
#' ### check that x[2,2,2] mean is correct
#' check = x[1:3,1:3,1:3]
#' ## masking
#' vals = check[abs(check) < 1]
#' m = mean(vals)
#' all.equal(m, mean.x[2,2,2])
#' n = length(vals)
#' v = var(vals) * (n-1)/n
#' var.x[2,2,2]
#' all.equal(v, var.x[2,2,2])
#'
local_moment <- function(
  image,
  window = NULL,
  nvoxels=NULL,
  moment,
  mask = NULL,
  only.mask = FALSE,
  center=is.null(mean_image),
  invert = FALSE,
  # the mean
  mean_image=NULL, # mean image to be subtracted.  If not supplied
  # local_moment_edge is run with mom = 1
  na.rm=TRUE, # remove NAs from the image (mask set to 0),
  remask = TRUE,  # set areas outside of mask to 0
  ...
  ) {

  is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) {
    abs(x - round(x)) < tol
  }

  ## if mask is null - whole image
  dimg = dim(image)[1:3]
  if (is.null(mask)) {
    mask = array(1, dim=dimg)
  }

  lmom = length(moment)
  lcent = length(center)
  stopifnot(lmom == lcent)

  neigh = get_neighbors(image, window = window, nvoxels=nvoxels,
                mask = mask, ...)

  neighbors = neigh$neighbors
  indices = neigh$indices

  mn = rowMeans(neighbors, na.rm=na.rm)
  img_list = vector("list", length = lmom)
  for (imom in seq(lmom)){
    cent = center[imom]
    mom = moment[imom]
    moment_image = neighbors
    if (isTRUE(cent)){
      if (mom != 1) moment_image = moment_image - mn
    }

    moment_image = moment_image^mom
    moment_image = rowMeans(moment_image, na.rm=na.rm)

    realpow <- function(x, pow) {
      sgn = sign(x)
      x = abs(x)
      x = x^{pow}
      x = sgn * x
    }
    ### put on same scale
    if (invert) {
      moment_image = realpow(moment_image, pow = 1/mom)
    }

    moment_image = array(moment_image, dim=dimg)
    if (remask) {
      moment_image = moment_image * mask
    }
    if ( inherits(image, "nifti") ){
      moment_image = niftiarr(image, moment_image)
      moment_image = datatyper(moment_image,
                              datatype= convert.datatype()$FLOAT32,
                              bitpix= convert.bitpix()$FLOAT32)
    }
    img_list[[imom]] = moment_image
  }

  # moment <- array(moment, dim = dim(image))
  return(img_list)
}


#' @title Calculate Moments of Neighborhood
#'
#' @description This function calculates Local Moments (mean, standard deviation, skew) for an array.
#' @param image input image
#' @param window window (in width) for the neighborhood
#' @param nvoxels window (in voxels) for the neighborhood 1 results in a 3x3 cube
#' @param mask array or object of class nifti of same size as image
#' @param rm.value remove the voxel itself before taking the moment
#' @param check.wrap Logical - check wrapround and put \code{rep.value} in
#' for the wrapped values
#' @param rep.value Replace wrapped values (edge of image) to this value
#' @param ... Not used
#' @importFrom magic ashift
#' @export
#' @return Array with same dimensions as image
#' @examples
#' x = array(rnorm(1000), dim=c(10, 10, 10))
#' neigh = get_neighbors(x, nvoxels = 1)
get_neighbors <- function(
  image,
  window = NULL,
  nvoxels=NULL,
  mask = NULL,
  rm.value = FALSE, # remove the voxel itself before returning
  check.wrap = TRUE,
  rep.value = 0,
  ...
) {

  is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) {
    abs(x - round(x)) < tol
  }

  ## if mask is null - whole image
  dimg = dim(image)[1:3]
  if (is.null(mask)) {
    mask = array(1, dim=dimg)
  }

  mypermutations = function(winds){
    windlist = lapply(1:3, function(x) winds)
    eg = expand.grid(windlist)
    eg = eg[ order(eg$Var1, eg$Var2, eg$Var3), ]
    eg = as.matrix(eg)
    rownames(eg) = NULL
    colnames(eg) = NULL
    eg
  }

  ### initalize array
  ### different ways of parameterizing the "window"
  if (is.null(nvoxels)) {
    stopifnot(!is.null(window))
    stopifnot(length(window) == 1)
    if ((window %% 2) != 1) {
      stop("window must be odd number")
    }
    winds = (-window/2 + .5):(window/2 - .5)
#     indices <- gtools::permutations(window, 3,
#                             v= (-window/2 + .5):(window/2 - .5),
#                             repeats.allowed=TRUE)
  }  else {
    stopifnot(is.wholenumber(nvoxels))
    stopifnot(is.null(window))

    winds <- (-nvoxels):nvoxels
#     indices <- gtools::permutations(length(winds), 3, v = winds,
#                             repeats.allowed=TRUE)
  }
  indices = mypermutations(winds)

  if (rm.value){
    allzero = apply(indices == 0, 1, all)
    indices = indices[!allzero,]
  }

  image = image * mask

  nruns <- nrow(indices)
  mat = matrix(NA, nrow=prod(dimg), ncol=nruns)
  for ( i in 1:nruns){
    shifter <- ashift(image, indices[i,])
    mat[, i] = shifter
  }

  if (check.wrap){
    tall.ind = t(as.matrix(expand.grid(dim1=seq(dimg[1]),
                                      dim2=seq(dimg[2]),
                                      dim3=seq(dimg[3]))))
    dimg.mat = matrix(dimg, ncol=prod(dimg), nrow = 3)
    for ( i in 1:nruns){
      all.ind2 = tall.ind + indices[i,]
      outside = {all.ind2 < 1} + all.ind2 > dimg.mat
      any.outside = colSums(outside) > 0
      mat[any.outside,i] = rep.value
    }
  }

  return(list(neighbors=mat, indices = indices))
}


#' @title Calculate Mean Image by FFT
#'
#' @description This function calculates the mean image using fft
#' @param x 3D array
#' @param nvoxels window (in voxels) for the neighborhood
#' (1 results in a 3x3 cube)
#' @param shift Should results be shifted back?
#' @param verbose Diagnostic outputing
#' @importFrom magic ashift
#' @importFrom stats fft
#' @export
#' @return Array of size of x
mean_image = function(x, nvoxels, shift = TRUE, verbose = TRUE){

  is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) {
    abs(x - round(x)) < tol
  }
  stopifnot(length(nvoxels) %in% c(1, 3))

  if (length(nvoxels) == 1){
    stopifnot(is.wholenumber(nvoxels))
    n <- length((-nvoxels):nvoxels)

    h = array(1, dim=c(n, n, n))
    h = h / length(h)
  }
  if (length(nvoxels) == 3){
    dims = (nvoxels*2) + 1
    stopifnot(is.wholenumber(dims))
    h = array(1, dim=dims)
    h = h / length(h)
  }
  if (verbose){
    cat(paste0("Dimension of kernel is ", paste(dim(h), collapse = "x"), "\n"))
  }

  # % x - 3dim matrix
  # % h - smoothing kernel - 3dim matrix size(h) =< size(x)

  x[is.na(x)]=0;

  dx = dim(x)
  RX = dx[1];
  CX = dx[2];
  SX = dx[3];

  dh = dim(h)
  RH = dh[1];
  CH = dh[2];
  SH = dh[3];

  z = array(0, dim = dx);
  z[1:RH,1:CH,1:SH] = h;

  if (verbose){
    cat("Creating Kernel fft\n")
  }
  H = fft(z);
  rm(z)
  if (verbose){
    cat("Creating Image fft\n")
  }
  X = fft(x);
  rm(x)
  if (verbose){
    cat("Convolution\n")
  }
  y = H * X
  rm(H)
  rm(X)
  y = fft(y, inverse = TRUE)
  y = y / length(y)
  y = Re(y)

#   print(ceiling(RH/2))
#   print(ceiling(CH/2))
#   print(ceiling(SH/2))
  #   [p+1:m 1:p]
  if (verbose){
    cat("Shifting\n")
  }
  if (shift) {
    y = ashift(y, v = -(ceiling(dim(h)/2)-1))
  }

  #   y=circshift(y,c(-ceiling(RH/2), -ceiling(CH/2), -ceiling(SH/2)))
  # compensation for group delay
  return(y)
}