#' @title Predict ICH Images

#' @description This function will take the data.frame of predictors and

#' predict the ICH voxels from the model chosen.

#'

#' @param df \code{\link{data.frame}} of predictors.  If \code{multiplier}

#' column does not exist, then \code{\link{ich_candidate_voxels}} will

#' be called

#' @param nim object of class \code{\link{nifti}}, from

#' \code{\link{make_predictors}}

#' @param model model to use for prediction,

#' either the random forest (rf) or logistic

#' @param verbose Print diagnostic output

#' @param native Should native-space predictions be given?

#' @param native_img object of class \code{\link{nifti}}, which

#' is the dimensions of the native image

#' @param transformlist Transforms list for the transformations back to native space.

#' NOTE: these will be inverted.

#' @param interpolator Interpolator for the transformation back to native space

#' @param native_thresh Threshold for re-thresholding binary mask after

#' interpolation

#' @param shiny Should shiny progress be called?

#' @param model_list list of model objects, used mainly for retraining

#' but only expert use.

#' @param smoothed_cutoffs A list with an element

#' \code{mod.dice.coef}, only expert use.

#' @param outfile filename for output file.

#' We write the smoothed, thresholded image.  If \code{native = TRUE},

#' then the file will be native space, otherwise in registered

#' space

#' @param ... Additional options passsed to \code{\link{ich_preprocess}}

#'

#' @return List of output registered and native space

#' prediction/probability images

#' @importFrom neurobase remake_img

#' @importFrom extrantsr ants_bwlabel

#' @import randomForest

#' @seealso \code{\link{ich_candidate_voxels}}

#' @export

ich_predict = function(df,

                       nim,

                       model = c("rf", "logistic", "big_rf"),

                       verbose = TRUE,

                       native = TRUE,

                       native_img = NULL,

                       transformlist = NULL,

                       interpolator = NULL,

                       native_thresh = 0.5,

                       shiny = FALSE,

                       model_list = NULL,

                       smoothed_cutoffs = NULL,

                       outfile = NULL,

                       ...) {

  # if (!have_matlab()) {

  #   stop("MATLAB Path not defined!")

  # }

  if (is.null(outfile)) {

    outfile = tempfile(fileext = ".nii.gz")

  }

  cn = colnames(df)

  if (!("multiplier" %in% cn)) {

    df$multiplier = ich_candidate_voxels(df)

  }

  df$Y = NULL

  cc = complete.cases(df)

  if (!all(cc)) {

    warning("NAs or missing in DF, removing")

    for (icn in seq(ncol(df))) {

      x = df[, icn]

      if (!(class(x) %in% c("factor", "character"))) {

        x[ !is.finite(x) ] = 0

      }

      df[, icn] = x

    }

  }

  msg = "# Making Prediction"

  if (verbose) {

    message(msg)

  }

  if (shiny) {

    shiny::incProgress(message = msg)

  }

  env = as.environment("package:ichseg")

  # Getting modlist for model and cutoff

  if (is.null(model_list)) {

    modlist.name = paste0(model, "_modlist")

    modlist = env[[modlist.name]]

  } else {

    modlist = model_list

  }

  mod = modlist$mod

  cutoff = modlist$mod.dice.coef[1, "cutoff"]

  rm(list = c("modlist"))

  # Getting smoothed cutoff

  if (is.null(smoothed_cutoffs)) {

    smoothed_name = paste0("smoothed_", model, "_cutoffs")

    scutoffs = env[[smoothed_name]]

  } else {

    scutoffs = smoothed_cutoffs

  }

  smoothed_cutoff = scutoffs$mod.dice.coef[1, "cutoff"]

  rm(list = c("scutoffs", "smoothed_name"))

  p = switch(model,

             rf = predict(mod,

                          newdata = df[ df$multiplier, ],

                          type = "prob")[, "1"],

             big_rf = predict(mod,

                              newdata = df[ df$multiplier, ],

                              type = "prob")[, "1"],

             logistic = predict(mod,

                                df[ df$multiplier, ],

                                type = "response"))

  msg = "# Making Prediction Image"

  if (verbose) {

    message(msg)

  }

  nim = check_nifti(nim)

  mult_img = niftiarr(nim, df$multiplier)

  # p = predict(mod, df[ df$multiplier, ], type = "response")

  pimg = remake_img(p,

                    nim,

                    mult_img)

  mask = niftiarr(nim, df$mask)

  pimg = mask_img(pimg, mask)

  msg = "# Smoothing Image"

  if (verbose) {

    message(msg)

  }

  sm.pimg  = mean_image(pimg,

                        nvoxels = 1,

                        verbose = verbose)

  sm.pimg[abs(sm.pimg) <

            .Machine$double.eps ^ 0.5 ] = 0

  sm.pimg = niftiarr(nim, sm.pimg)

  sm.pimg[is.na(sm.pimg)] = 0

  sm.pred = sm.pimg > smoothed_cutoff

  pred = pimg > cutoff

  msg = "# Connected Components"

  if (verbose) {

    message(msg)

  }

  # cc = spm_bwlabel(pred, k = 100)

  # scc = spm_bwlabel(sm.pred, k = 100)

  cc = ants_bwlabel(img = pred, k = 100, binary = TRUE)

  scc = ants_bwlabel(img = sm.pred, k = 100, binary = TRUE)

  ##############################################################

  # Back to Native Space!

  ##############################################################

  res = list(

    prediction_image = cc,

    smoothed_prediction_image = scc,

    probability_image = pimg,

    smoothed_probability_image = sm.pimg)

  ##############################################################

  # Inverted!

  ##############################################################

  native_res = NULL

  if (native) {

    msg = "# Projecting back to Native Space"

    if (verbose) {

      message(msg)

    }

    stopifnot(!is.null(interpolator))

    stopifnot(!is.null(transformlist))

    native_res = lapply(res, function(x){

      ants_apply_transforms(fixed = native_img,

                            moving = x,

                            transformlist = transformlist,

                            interpolator = interpolator,

                            whichtoinvert = c(1)

      )

    })

    native_res$smoothed_prediction_image = neurobase::datatyper(

      native_res$smoothed_prediction_image > native_thresh

    )

    native_res$prediction_image = neurobase::datatyper(

      native_res$prediction_image > native_thresh

    )

    writenii(native_res$smoothed_prediction_image, outfile)

  } else {

    writenii(res$smoothed_prediction_image, outfile)

  }

  res$cutoff = cutoff

  res$smoothed_cutoff = smoothed_cutoff

  L = list(registered_prediction = res)

  L$native_prediction = native_res

  L$outfile = outfile

  return(L)

}