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

}