#' @title Coarse Masking of Image from Z-slicding
#' @description Calculates an Ellipsoid at each slice then stacks them
#'
#' @param x Binary \code{antsImage}
#'
#' @return Array or \code{antsImage} object
#' @export
#' @importFrom ANTsRCore is.antsImage
#' @importFrom stats cov
ellipsoid_body = function(x) {
  bb = as.array(x)
  dimg = dim(bb)
  number_slices = dimg[3]
  dslice = dimg[1:2]
  omat = array(FALSE, dim = dimg)

  # get values that are true for ellipse
  ind = which(bb > 0, arr.ind = TRUE)
  L = lapply(dslice, seq)
  # same for all slices
  X = as.matrix(expand.grid(L))

  # need df for split
  ind = as.data.frame(ind)
  ps = split(ind, ind[,3])
  ps = lapply(ps, function(x){
    as.matrix(x[, 1:2])
  })
  # Getting mu/sigma for ellipse
  vals = lapply(ps, function(x){
    L = list(mu = colMeans(x),
             Sigma = cov(x))
  })
  # Taken from SIBER, but faster because no looping over solve
  pointsToEllipsoid2 = function(X, Sigma, mu) {
    eig <- eigen(Sigma)
    SigSqrt = eig$vectors %*% diag(sqrt(eig$values)) %*% t(eig$vectors)
    xtx = t(SigSqrt) %*% SigSqrt
    ss = solve(xtx)
    xx = t(t(X) - mu)
    Z = ss %*% t(SigSqrt) %*% t(xx)
    Z = t(Z)
  }
  # Taken from SIBER
  ellipseInOut2 = function(Z, p = 0.95, radius = NULL) {
    if (is.null(radius)) {
      radius <- stats::qchisq(p, df = ncol(Z))
    }
    inside <- rowSums(Z^2) < radius
    return(inside)
  }
  for ( i in seq(number_slices)) {
    x = vals[[i]]
    Z = pointsToEllipsoid2(X, Sigma = x$Sigma, mu = x$mu)
    inside = ellipseInOut2(Z = Z, p = 0.95)
    mat = array(inside, dim = dslice)
    omat[,,i] = mat
  }
  if (is.antsImage(x)) {
    omat = as.antsImage(omat, reference = x)
  }
  return(omat)
}