#' @title Make CT Predictors
#'
#' @description Create a set of predictors for ICH segmentation
#' for CT
#' @param img Filename of image intensities
#' @param mask Filename of brain mask
#' @param roi Filename of ROI for Y
#' @param nvoxels Voxel neighborhood
#' @param moments Moments of neighborhood to take
#' @param center Center the moments
#' @param lthresh Lower threshold for neighborhood setting
#' @param uthresh Upper threshold for neighborhood setting
#' @param sigmas Sigma values for Gaussian smoothed images (in mm)
#' @param save_imgs Logical to save all images that are created as
#' predictors
#' @param outdir Output directory of saved images, needs to be set
#' if \code{save_imgs = TRUE}
#' @param stub Basename to write image names if \code{save_imgs = TRUE}
#' @param overwrite If \code{save_imgs} is \code{TRUE},
#' then should
#' the files be overwritten?  If not, then files will be read
#' instead instead of code being re-run.
#' @param template.file Template to register to (CT Template)
#' @param mean.img Mean image in template space for z-scoring
#' @param sd.img SD image in template space for z-scoring
#' @param zscore.typeofTransform type of transform for z-scoring
#' @param zscore.interpolator type of interpolator for z-scoring
#' @param flip.typeofTransform type of transform for flipped difference
#' @param flip.interpolator type of interpolator for flipped difference
#' @param low_thresh Threshold for forcing values to zero
#' @param verbose Logical indicator if output messages should be
#' printed
#' @param shiny Should shiny progress be called?
#' @param erode_mask Should the brain mask be eroded?
#' @param ... options passed to \code{\link{get_neighbors}}
#' @importFrom neurobase readnii checkimg zscore_img finite_img
#' @importFrom fslr fslerode fsl_smooth
#' @importFrom oro.nifti zero_trans cal_img voxdim pixdim convert.datatype convert.bitpix
#' @importFrom extrantsr zscore_template otropos reg_flip perona_malik
#' @importFrom stats sd quantile predict complete.cases
#' @export
#' @return List of a data.frame of Predictors and set of
#' indices to
#' keep in mask and an empty nifti object for plotting.
#' Also the number of voxels of the roi that were not in the
#' mask
make_predictors <- function(img, mask, roi = NULL,
                            nvoxels = 1, moments = 1:4,
                            center = c(FALSE, TRUE, TRUE, TRUE),
                            lthresh = 40, uthresh = 80,
                            sigmas = c(5, 10, 20),
                            save_imgs = TRUE,
                            outdir = NULL,
                            stub = NULL,
                            overwrite = FALSE,
                            template.file = system.file(
                              "scct_unsmooth_SS_0.01.nii.gz",
                              package = "ichseg"),
                            mean.img = system.file(
                              "Mean_Image.nii.gz",
                              package = "ichseg"),
                            sd.img = system.file(
                              "SD_Image.nii.gz",
                              package = "ichseg"),
                            zscore.typeofTransform = "SyN",
                            zscore.interpolator = "Linear",
                            flip.typeofTransform = "Affine",
                            flip.interpolator = "LanczosWindowedSinc",
                            low_thresh = 1e-13,
                            verbose= TRUE,
                            shiny = FALSE,
                            erode_mask = TRUE,
                            ...) {
  make_fname = function(addstub){
    fname = addstub
    fname = file.path(outdir, paste0(stub, "_", fname, ".nii.gz"))
    fname
  }

  write_img = function(arr, addstub){
    fname = addstub
    fname = file.path(outdir, paste0(stub, "_", fname))
    if ( !inherits( arr, "nifti")){
      mom = nim
      mom@.Data = array(arr, dim = dim(mom))
    } else {
      mom = arr
    }
    mom = cal_img(mom)
    mom = zero_trans(mom)
    mom = datatyper(mom,
                    datatype = convert.datatype()$FLOAT32,
                    bitpix = convert.bitpix()$FLOAT32)
    writenii(mom, filename = fname)
  }


  make_moment_list = function(){
    mom.exist = sapply(moments, function(moment){
      addstub = paste0("moment", moment)
      fname = make_fname(addstub)
      file.exists(fname)
    })

    ### if all data does not exist or rewrite - remake data
    if (!all(mom.exist) | overwrite) {
      msg = "# Running Local_Moment"
      if (verbose){
        message(msg)
      }
      if (shiny) {
        shiny::incProgress(message = msg, amount = 0.02)
      }
      moms = local_moment(
        img,
        mask = mask,
        nvoxels = nvoxels,
        moment = moments, center = center,
        invert = FALSE,
        ...)
      if (save_imgs) {
        for (imom in seq_along(moments)) {
          moment = moments[imom]
          addstub = paste0("moment", moment)
          fname = make_fname(addstub)
          mom.img = moms[[imom]]
          write_img(mom.img, addstub)
        }
      }
    }
    ### if all data exists and not to rewrite - just read in
    if (all(mom.exist) & !overwrite) {
      moms = lapply(moments, function(moment){
        addstub = paste0("moment", moment)
        fname = make_fname(addstub)
        mom.img = readnii(fname, reorient = FALSE)
        return(mom.img)
      })
    }
    msg = "# Creating Moment Matrix"
    if (verbose) {
      message(msg)
    }
    if (shiny) {
      shiny::incProgress(message = msg, amount = 0.02)
    }
    mat = matrix(NA,
                 ncol = length(moms),
                 nrow = prod(dim(mask))
    )
    for (imom in seq_along(moms)) {
      mat[, imom] = c(moms[[imom]])
    }
    # moms = sapply(moms, c)
    rm(list = "moms"); gc(); gc();
    colnames(mat) = paste0("moment", moments)
    mat = as.data.frame(mat)
    return(mat)
  }

  img.fname = checkimg(img)
  img = check_nifti(img)

  mask.fname = checkimg(mask)
  orig.mask = mask = check_nifti(mask)

  if (!is.null(roi)) {
    #     roi.fname = roi
    roi = check_nifti(roi)
  }

  # stopifnot(class(roi) == "character")
  #   stopifnot(class(img) == "character")
  #   stopifnot(class(mask) == "character")


  if (save_imgs){
    stopifnot(!is.null(outdir))
    stopifnot(!is.null(stub))
  }
  if (is.null(outdir)){
    outdir = tempdir()
  }
  if (is.null(stub)){
    stub = nii.stub(img.fname, bn = TRUE)
  }

  msg = "# Reading Images"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  # img = readnii(img.fname, reorient= FALSE)
  orig.img = img
  # stub = nii.stub(img.fname, bn=TRUE)

  nim = niftiarr(img, array(NA, dim = dim(orig.img)))
  nim = datatyper(nim,
                  datatype = convert.datatype()$FLOAT32,
                  bitpix = convert.bitpix()$FLOAT32)
  dimg = dim(nim)
  vdim = voxdim(nim)

  #   if (is.null(roi)){
  #     roi = niftiarr(img, array(0, dim = dim(nim)))
  #   }

  msg = "# Eroding Mask"
  addstub = "usemask"
  fname = make_fname(addstub)

  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }

  if (file.exists(fname) & !overwrite){
    mask = readnii(fname, reorient = FALSE)
  } else {
    if (erode_mask) {
      # erode the mask
      mask = fslerode(file = mask.fname,
                      kopts = "-kernel box 3x3x1",
                      reorient = FALSE, retimg = TRUE,
                      verbose = verbose > 1)
    } else {
      mask = readnii(mask.fname)
    }
    #### may add this - think about it
    #     mask = fslfill(mask, bin=TRUE, retimg = TRUE)
    mask = mask > 0
    mask = datatyper(mask)
    if (save_imgs){
      write_img(mask, addstub)
    }
  }

  if (sum(mask) == 0) {
    msg = paste0(
      "Eroded mask is empty! Something went wrong with eroding ",
      "or Skull stripping")
    stop(msg)
  }
  orig.masked.img = mask_img(orig.img, orig.mask)
  masked.img = mask_img(orig.img, mask)
  if (!is.null(roi)){
    miss.roi = sum(mask == 0 & roi > 0 )
  } else {
    miss.roi = NULL
  }


  keep.ind = which(mask > 0)

  msg = "# Getting Moments"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  ################################################
  # Making Moment Images
  ################################################

  df = make_moment_list()
  ###########
  # Creating Skew/Kurtosis
  ###########
  df$skew = df$moment3/df$moment2 ^ {3/2}
  df$kurtosis = df$moment4/df$moment2 ^ {2}
  df$moment3 = df$moment4 = NULL

  ###########
  # Writing Skew
  ###########
  addstub = "skew"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite) {

  } else {
    skew = niftiarr(img, df$skew)
    if (save_imgs) {
      write_img(skew, addstub)
    }
    rm(list = c("skew")); gc(); gc();
  }

  ###########
  # Writing Kurtosis
  ###########
  addstub = "kurtosis"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite) {

  } else {
    kurtosis = niftiarr(img, df$kurtosis)
    if (save_imgs) {
      write_img(kurtosis, addstub)
    }
    rm(list = c("kurtosis")); gc(); gc();
  }

  ################################################
  # Making Percent threshold image
  ################################################
  msg = paste0("# Getting thresholded from ",
               lthresh, " to ", uthresh, "\n")
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = paste0("thresh_", lthresh, "_", uthresh)
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    thresh_img = readnii(fname, reorient=FALSE)
  } else {
    thresh_img = niftiarr(img, img >= lthresh & img <= uthresh)
    if (save_imgs){
      write_img(thresh_img, addstub)
    }
  }

  df$value = c(masked.img)
  df$thresh = c(thresh_img)

  cn = colnames(df)

  ################################################
  # Making Z-score Images
  ################################################
  msg = "# Getting Z-scored images"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  zimgs = lapply(1:3, function(i){
    addstub = paste0("zscore", i)
    fname = make_fname(addstub)
    if (file.exists(fname) & !overwrite){
      img = readnii(fname, reorient=FALSE)
    } else {
      img = zscore_img(masked.img, mask = mask, margin = i)
      if (save_imgs){
        write_img(img, addstub)
      }
    }
    return(img)
  })
  zimgs = sapply(zimgs, c)
  colnames(zimgs) = paste0("zscore", 1:3)
  zimgs = as.data.frame(zimgs)
  for (i in 1:3) {
    df[,  paste0("zscore", i)] = zimgs[, paste0("zscore", i)]
  }
  rm(list = c("zimgs")); gc(); gc();

  #   wmean2 = function(img, mask, trim = 0.2){
  #     x = img[ mask == 1]
  #     mn = psych::winsor.mean(x, trim = trim)
  #     s = psych::winsor.sd(x, trim = trim)
  #     z = (x-mn)/s
  #     img[mask == 1] = z
  #     img[mask == 0] = 0
  #     img = cal_img(img)
  #     return(img)
  #   }

  wmean = function(img, mask, trim = 0.2){
    x = img[ mask == 1 ]
    stopifnot(length(trim) == 1)
    stopifnot(trim > 0)
    stopifnot(trim <= 0.5)
    qtrim <- quantile(x,
                      c(trim, 0.5, 1 - trim),
                      na.rm = TRUE)
    xbot <- qtrim[1]
    xtop <- qtrim[3]

    if (trim < 0.5) {
      x[x < xbot] <- xbot
      x[x > xtop] <- xtop
    } else {
      x[!is.na(x)] <- qtrim[2]
    }

    mn = mean(x, na.rm=TRUE)
    s = sd(x, na.rm=TRUE)
    img = (img-mn)/s
    img = mask_img(img, mask)
    img = finite_img(img)
    return(img)
  }

  msg = "# Getting Winsorized Image"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = paste0("win_z")
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    wmean_img = readnii(fname, reorient=FALSE)
  } else {
    wmean_img = wmean(orig.img, mask = mask, trim = 0.2)
    if (save_imgs){
      write_img(wmean_img, addstub)
    }
  }
  df$win_z = c(wmean_img)
  rm(list = c("wmean_img")); gc(); gc();

  ################################################
  # Making Percent threshold image
  ################################################
  msg = paste0("# Getting Percent thresholded from ",
               lthresh, " to ", uthresh)
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = paste0("pct_thresh_", lthresh, "_", uthresh)
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    pct_thresh = readnii(fname, reorient=FALSE)
  } else {
    pct_thresh = local_moment(thresh_img, mask = mask,
                              nvoxels = nvoxels,
                              moment = 1, center = FALSE,
                              ...)[[1]]
    if (save_imgs){
      write_img(pct_thresh, addstub)
    }
  }

  df$pct_thresh = c(pct_thresh)
  rm(list = c("thresh_img", "pct_thresh")); gc(); gc();


  ################################################
  # Making Probability
  ################################################
  msg = "# Getting Top Probability Segmentation from Atropos"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = paste0("prob_img")
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    prob_img = readnii(fname, reorient = FALSE)
  } else {
    window.masked.img = window_img(masked.img, window = c(0, 100))
    seg = otropos( window.masked.img, i = "KMeans[4]", verbose = verbose > 1)
    prob_img = seg$probabilityimages[[3]] + seg$probabilityimages[[4]]
    rm(list = c("seg")); gc(); gc();
    if (save_imgs){
      write_img(prob_img, addstub)
    }
  }

  df$prob_img = c(prob_img)
  rm(list = c("prob_img")); gc(); gc();


  ################################################
  # Making Percent that are 0
  ################################################
  msg = "# Getting Percent 0"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  ###### changed to masked.img
  thresh_0 = niftiarr(masked.img, masked.img == 0)
  addstub = "pct_zero_neighbor"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    pct_zero_neighbor = readnii(fname, reorient=FALSE)
  } else {
    pct_zero_neighbor = local_moment(thresh_0, mask = NULL,
                                     nvoxels = nvoxels,
                                     moment = 1, center = FALSE,
                                     ...)[[1]]
    pct_zero_neighbor = mask_img(pct_zero_neighbor, mask)
    if (save_imgs){
      write_img(pct_zero_neighbor, addstub)
    }
  }
  df$pct_zero_neighbor = c(pct_zero_neighbor)
  rm(list = c("pct_zero_neighbor", "thresh_0")); gc(); gc();

  msg = "# Getting Any 0 Neighbors"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = "any_zero_neighbor"
  df$any_zero_neighbor = (df$pct_zero_neighbor > 0) *1

  ################################################
  # Making Distance to centroid
  ################################################
  msg = "# Getting Distance to centroid"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = "dist_centroid"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite){
    dist.img = readnii(fname, reorient=FALSE)
  } else {
    centroid = t(which(mask > 0, arr.ind=TRUE))
    all.ind = expand.grid(lapply(dimg, seq))
    colnames(all.ind) = paste0("dim", seq(length(dimg)))
    all.ind = t(as.matrix(all.ind))
    all.ind = all.ind * vdim
    centroid = centroid * vdim
    dist.img = t(all.ind - rowMeans(centroid))
    rm(list = c("all.ind")); gc(); gc();

    dist.img = sqrt(rowSums(dist.img^2))
    dist.img = niftiarr(img, array(dist.img, dim =dimg))
    dist.img = mask_img(dist.img, mask)
    dist.img = datatyper(dist.img,
                         datatype= convert.datatype()$FLOAT32,
                         bitpix= convert.bitpix()$FLOAT32)
    if (save_imgs){
      write_img(dist.img, addstub)
    }
  }
  df$dist_centroid = c(dist.img)
  rm(list = c("dist.img")); gc(); gc();


  ################################################
  # Making Distance to centroid
  ################################################
  msg = "# Perona Malik Smoother"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = "perona_malik"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite) {
    dist.img = readnii(fname, reorient = FALSE)
  } else {
    pm_img = extrantsr::perona_malik(
      masked.img, n_iter = 10,
      conductance = 5)
    if (save_imgs) {
      write_img(pm_img, addstub)
    }
    rm(list = c("seg")); gc(); gc();
  }
  df$perona_malik = c(pm_img)
  rm(list = c("pm_img")); gc(); gc();

  ################################################
  # Making 10mm and 20mm smoother
  ################################################
  make_smooth_img = function(sigma){
    msg = paste0("# Getting Smooth ", sigma)
    if (verbose) {
      message(msg)
    }
    if (shiny) {
      shiny::incProgress(message = msg, amount = 0.02)
    }
    addstub = paste0("smooth", sigma)
    if (save_imgs){
      fname = make_fname(addstub)
    } else {
      fname = tempfile()
      if (file.exists(fname)){
        file.remove(fname)
      }
    }
    if (file.exists(fname) & !overwrite) {
      smooth.img = readnii(fname, reorient = FALSE)
    } else {
      smooth.mask = fslr::fsl_smooth(
        file = mask,
        sigma = sigma,
        mask = NULL,
        smooth_mask = FALSE,
        verbose = verbose > 1)
      smooth.img = fslsmooth(img.fname, sigma=sigma,
                             mask = mask, retimg = TRUE,
                             outfile = fname,
                             verbose = verbose > 1)
    }
    return(c(smooth.img))
  }
  #   df$smooth2 = make_smooth_img(sigma=2)
  #   df$smooth5 = make_smooth_img(sigma=5)
  smooths = sapply(sigmas, make_smooth_img)
  colnames(smooths) = paste0("smooth", sigmas)
  df = cbind(df, smooths)
  rm(list = c("smooths")); gc(); gc();



  msg = "# Z-score to template"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = "zscore_template"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite) {
    zscore = readnii(fname, reorient = FALSE)
  } else {
    zscore = extrantsr::zscore_template(img = orig.masked.img,
                                        template.file = template.file,
                                        mean.img = mean.img,
                                        sd.img = sd.img,
                                        typeofTransform = zscore.typeofTransform,
                                        interpolator = zscore.interpolator,
                                        verbose = verbose > 1)
    if (save_imgs){
      write_img(zscore, addstub)
    }
  }
  df$zscore_template = c(zscore)
  rm(list = c("zscore")); gc(); gc();

  msg = "# Flipped Difference"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }
  addstub = "flipped_value"
  fname = make_fname(addstub)
  if (file.exists(fname) & !overwrite) {
    flipped_value = readnii(fname, reorient=FALSE)
  } else {
    flipper = extrantsr::reg_flip(
      t1 = orig.masked.img,
      mask = orig.mask,
      template.file = template.file,
      typeofTransform = flip.typeofTransform,
      interpolator = flip.interpolator,
      verbose = verbose > 1)
    flipper = flipper$t1
    ##########################
    # Take difference
    ##########################
    flipped_value = orig.masked.img - flipper
    rm(list = c("flipper")); gc(); gc();
    if (save_imgs) {
      write_img(flipped_value, addstub)
    }
  }
  df$flipped_value = c(flipped_value)
  rm(list = c("flipped_value")); gc(); gc();

  msg = "# Thresholding small values"
  if (verbose) {
    message(msg)
  }
  if (shiny) {
    shiny::incProgress(message = msg, amount = 0.02)
  }


  for (icn in seq(ncol(df))) {
    x = df[, icn]
    if (!(class(x) %in% c("factor", "character"))) {
      x[ !is.finite(x) ] = 0
    }
    df[, icn] = x
  }

  df = as.matrix(df)
  low = abs(df) < low_thresh
  df[ low ] = 0

  df = data.frame(df, stringsAsFactors = FALSE)
  if (!is.null(roi)) {
    df$Y = c(roi)
  } else {
    df$Y = NA
  }
  df$mask = c(mask)


  return(list(df = df, keep.ind = keep.ind, nim = nim,
              miss.roi = miss.roi))
}