#' @title Fill image holes with dilation then erosion
#' @description This function calls \code{mean_image} to dilate an image, then calls
#' it again to erode it.
#' @param file (character) filename of image to be filled
#' @param outfile (character) name of resultant filled file
#' @param nvoxels (integer) Number of voxels to smooth over, creates vxvxv box.
#' @param zeropad (logical) Perform \code{zero_pad} before running.
#' @param remove.ends (logical) Remove top and bottom dilation.
#' @param tol (double) Tolerance for thresholding after \code{\link{fft}}.
#' @param refill (logical) Run \code{\link{fslfill}} after dilation/erosion.
#' @param retimg (logical) return image of class nifti
#' @param reorient (logical) If retimg, should file be reoriented when read in?
#' Passed to \code{\link{readNIfTI}}.
#' @param intern (logical) pass to \code{\link{system}}
#' @param verbose (logical) print out command before running
#' @param ... additional arguments passed to \code{\link{readNIfTI}}.
#' @return character or logical depending on intern
#' @importFrom neurobase zero_pad
#' @note This function binarizes the image before running.
#' @export
dil_ero = function(file,
                   outfile = NULL,
                   nvoxels = 3,
                   zeropad = TRUE,
                   remove.ends = FALSE,
                   tol = .Machine$double.eps^0.5,
                   refill = TRUE,
                   retimg = FALSE,
                   reorient = FALSE,
                   intern=TRUE, verbose = TRUE,
                   ...){
  have.outfile = TRUE

  if (retimg){
    if (is.null(outfile)) {
      outfile = tempfile()
      have.outfile = FALSE
    }
  } else {
    stopifnot(!is.null(outfile))
  }

  file = check_nifti(file, reorient = reorient)
  img = file > 0
  dimg = dim(file)
  ##### should make for all max
  ind = which(img > 0, arr.ind = TRUE)
  ind = ind[ (ind[, "dim3"] %in% c(1, dimg[3])) |
               (ind[, "dim1"] %in% c(1, dimg[1])) |
               (ind[, "dim2"] %in% c(1, dimg[2])) ,]
  nind = nrow(ind)

  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)
    kdim = c(n, n, n)
  }
  if (length(nvoxels) == 3){
    kdim = (nvoxels*2) + 1
    stopifnot(is.wholenumber(kdim))
  }

  if (zeropad) {
    img = zero_pad(img, kdim = kdim)
  }
  dil = mean_image(img, nvoxels = nvoxels, verbose = verbose) > tol
  rm(img)
  dil = mean_image(1 - dil, nvoxels = nvoxels, verbose = verbose)
  dil = dil > tol
  dil = 1 - dil
  if (zeropad) {
    dil = zero_pad(dil, kdim = kdim, invert = TRUE)
  }
  dil = niftiarr(file, dil)
  rm(file)

  if (remove.ends) {
    #### making the ends correct - boundary problem
    dil@.Data[,,1] = array(0, dim=dimg[c(1,2)])
    dil@.Data[,,dimg[3]] = array(0, dim=dimg[c(1,2)])

    dil@.Data[,1,] = array(0, dim=dimg[c(1,3)])
    dil@.Data[,dimg[2],] = array(0, dim=dimg[c(1,3)])

    dil@.Data[1,,] = array(0, dim=dimg[c(2,3)])
    dil@.Data[dimg[1],,] = array(0, dim=dimg[c(2,3)])

    if (nind >0 ){
      dil@.Data[ ind ] = 1
    }
  }
  if (refill) {
    dil = fslfill(file = dil,
                  retimg=TRUE,
                  intern=intern, verbose = verbose)
  }
  dil = cal_img(dil)
  if (have.outfile){
    gzipped = grepl("gz$", get.imgext())
    writenii(dil, filename = outfile, gzipped = gzipped)
  }
  if (retimg){
    return(dil)
  }
  return(outfile)
}