#' @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)

}