#' Calculate first order radiomic features on a 2D or 3D array
#'
#' @param data Any 2D or 3D image (as matrix or array) to calculate first-order features
#' @param features = first level radiomic features to calculate
#'
#' @return Values from selected features
#' @importFrom stats quantile sd
#' @export
radiomics_first <- function(data,
                       features = c('mean', 'sd', 'skew', 'kurtosis', 'min', 'q1', 'median', 'q3', 'max','energy', 'rms', 'uniformity', 'entropy')){

  # Clean up data
  data <- as.vector(data)
  data <- data[!is.na(data)]

  # Moments
  if('mean' %in% features){
    mean_value <- mean(data)
  }else(mean_value = NULL)
  if('sd' %in% features){
    sd_value <- sd(data)
  }else(sd_value = NULL)
  if('skew' %in% features){
    skew_value <- skew(data)
  }else(skew_value = NULL)
  if('kurtosis' %in% features){
    kurtosis_value <- kurtosis(data)
  }else(kurtosis_value = NULL)


  # Quartiles
  if('min' %in% features){
    min_value <- min(data)
  }else(min_value = NULL)
  if('q1' %in% features){
    q1_value <- unname(quantile(data, probs = 0.25))
  }else(q1_value = NULL)
  if('median' %in% features){
    median_value <- median(data)
  }else(median_value = NULL)
  if('q3' %in% features){
    q3_value <- unname(quantile(data, probs = 0.75))
  }else(q3_value = NULL)
  if('max' %in% features){
    max_value <- max(data)
  }else(max_value = NULL)


  # Others
  if('energy' %in% features){
    energy_value <- energy(data)
  }else(energy_value = NULL)
  if('rms' %in% features){
    rms_value <- rms(data)
  }else(rms_value = NULL)
  if('uniformity' %in% features){
    uniformity_value <- uniformity(data)
  }else(uniformity_value = NULL)
  if('entropy' %in% features){
    entropy_value <- entropy(data)
  }else(entropy_value = NULL)


  featuresList <- list(
    mean = mean_value,
    sd = sd_value,
    skew = skew_value,
    kurtosis = kurtosis_value,
    min = min_value,
    q1 = q1_value,
    median = median_value,
    q3 = q3_value,
    max = max_value,
    energy = energy_value,
    rms = rms_value,
    uniformity = uniformity_value,
    entropy = entropy_value
  )
  if(length(features)==1){
    featuresList = unlist(featuresList[features])
  }

  return(featuresList)

}



skew <- function(data) {
  avg <- mean(data)
  SD  <- stats::sd(data)
  output <- mean(((data-avg)^3))/(SD)^3
  return(output)
}

kurtosis <- function(data) {
  avg <- mean(data)
  SD  <- stats::sd(data)
  output <- mean(((data-avg)^4))/(SD)^4
  return(output-3)
}

energy <- function(data) {
  output <- sum(data^2)
  return(output)
}

rms <- function(data) {
  output <- sqrt(sum(data^2)/length(data))
  return(output)
}

uniformity <- function(data) {
  output <- sum((table(data)/length(data))^2)
  return(output)
}

entropy <- function (data, base = 2)
{
  p <- table(data)/length(data)
  l <- logb(p, base)
  H <- sum(p * l)*-1
  return(H)
}