#' Calculate spatial radiomic features on a 2D or 3D array
#'
#' @param data Any 2D or 3D image (as matrix or array) to calculate spatial features
#' @param features = spatial radiomic features to calculate
#' @return Values from selected features
#' @importFrom abind abind
#' @export
radiomics_spatial <- function(data,
                              features = c('mi', 'gc', 'fd')){


  # Figure data dimension
  dimData <- length(dim(data))

  #2D
  if(dimData == 2){
    if('mi' %in% features){
      mi_value <- moran2D(data)
    }else(mi_value = NULL)
    if('gc' %in% features){
      gc_value <- geary2D(data)
    }else(gc_value = NULL)
    if('fd' %in% features){
      fd_value <- fd2D(data)
    }else(fd_value = NULL)

  }

  # 3D
  if(dimData == 3){

    if('mi' %in% features){
      mi_value <- moran3D(data)
    }else(mi_value = NULL)
    if('gc' %in% features){
      gc_value <- geary3D(data)
    }else(gc_value = NULL)
    if('fd' %in% features){
      fd_value <- NULL #fd3D(data) # still need to work on this
    }else(fd_value = NULL)

  }

  featuresList <- list(
    mi = mi_value,
    gc = gc_value,
    fd = fd_value
  )
  if(length(features)==1){
    featuresList = unlist(featuresList[features])
  }

  return(featuresList)

}


# Calculate Moran's I in 3D
moran3D <- function(mat){

  # Set up
  n <- length(mat[is.na(mat)==F])
  xbar <- mean(mat, na.rm=T)
  diff <- mat - xbar
  diff2 <- diff**2
  sum_diff2 <- sum(diff2, na.rm = T)

  # Dimension of matrix
  nx <- dim(diff)[1]
  ny <- dim(diff)[2]
  nz <- dim(diff)[3]

  # Rook shift
  left <- diff * abind(diff[-1,,], matrix(NA, nrow = ny, ncol = nz), along = 1)
  right <- diff * abind(matrix(NA, nrow = ny, ncol = nz), diff[-nx,,], along = 1)
  front <- diff * abind(diff[,-1,], matrix(NA, nrow = nx, ncol = nz), along = 2)
  back <- diff * abind(matrix(NA, nrow = nx, ncol = nz), diff[,-ny,], along = 2)
  up <- diff * abind(diff[,,-1], matrix(NA, nrow = nx, ncol = ny), along = 3)
  down <- diff * abind(matrix(NA, nrow = nx, ncol = ny), diff[,,-nz], along = 3)

  # Calcultate weights
  weights <- sum(length(left[is.na(left) == F])) +
    sum(length(right[is.na(right) == F])) +
    sum(length(front[is.na(front) == F])) +
    sum(length(back[is.na(back) == F])) +
    sum(length(up[is.na(up) == F])) +
    sum(length(down[is.na(down) == F]))

  # Change all "NAs" to zero, so we can efficiently sum matrices
  left[is.na(left)] <- 0
  right[is.na(right)] <- 0
  front[is.na(front)] <- 0
  back[is.na(back)] <- 0
  up[is.na(up)] <- 0
  down[is.na(down)] <- 0

  # Put info together to calculate MI
  sumv <- left + right + front + back + up + down
  local_mi <- sumv / sum_diff2 * n / weights
  global_mi <- sum(local_mi)

  return(global_mi)
}




# Calculate Moran's I in 2D
moran2D<-function(mat){

  # Set up
  n <- length(mat[is.na(mat)==F])
  xbar <- mean(mat, na.rm=T)
  diff <- mat - xbar
  diff2 <- diff**2
  sum_diff2 <- sum(diff2, na.rm = T)

  # Dimension of matrix
  nrow <- dim(diff)[1]
  ncol <- dim(diff)[2]

  # Rook shift
  up <- diff * rbind(diff[-1,],rep(NA,ncol))
  down <- diff * rbind(rep(NA,ncol),diff[-nrow,])
  left <- diff * cbind(diff[,-1],rep(NA,nrow))
  right <- diff * cbind(rep(NA,nrow),diff[,-ncol])

  # Calcultate weights
  weights <- sum(length(left[is.na(left) == F])) +
    sum(length(right[is.na(right) == F])) +
    sum(length(up[is.na(up) == F])) +
    sum(length(down[is.na(down) == F]))

  # Change all "NAs" to zero, so we can efficiently sum matrices
  left[is.na(left)] <- 0
  right[is.na(right)] <- 0
  up[is.na(up)] <- 0
  down[is.na(down)] <- 0

  # Put info together to calculate MI
  sumv <- left + right + up + down
  local_mi <- sumv / sum_diff2 * n / weights
  global_mi <- sum(local_mi)

  return(global_mi)
}




# Calculate Geary's C in 3D
geary3D<-function(mat){

  # Set up
  n <- length(mat[is.na(mat)==F])
  xbar <- mean(mat, na.rm=T)
  diff <- mat - xbar
  diff2 <- diff**2
  sum_diff2 <- sum(diff2, na.rm = T)

  # Dimension of matrix
  nx <- dim(diff)[1]
  ny <- dim(diff)[2]
  nz <- dim(diff)[3]

  # Rook shift
  left <- (mat - abind(mat[-1,,], matrix(NA, nrow = ny, ncol = nz), along = 1) )^2
  right <- (mat - abind(matrix(NA, nrow = ny, ncol = nz), mat[-nx,,], along = 1) )^2
  front <- (mat - abind(mat[,-1,], matrix(NA, nrow = nx, ncol = nz), along = 2) )^2
  back <- (mat - abind(matrix(NA, nrow = nx, ncol = nz), mat[,-ny,], along = 2) )^2
  up <- (mat - abind(mat[,,-1], matrix(NA, nrow = nx, ncol = ny), along = 3) )^2
  down <- (mat - abind(matrix(NA, nrow = nx, ncol = ny), mat[,,-nz], along = 3) )^2

  # Calcultate weights
  weights <- sum(length(left[is.na(left) == F])) +
    sum(length(right[is.na(right) == F])) +
    sum(length(front[is.na(front) == F])) +
    sum(length(back[is.na(back) == F])) +
    sum(length(up[is.na(up) == F])) +
    sum(length(down[is.na(down) == F]))

  # Change all "NAs" to zero, so we can efficiently sum matrices
  left[is.na(left)] <- 0
  right[is.na(right)] <- 0
  front[is.na(front)] <- 0
  back[is.na(back)] <- 0
  up[is.na(up)] <- 0
  down[is.na(down)] <- 0

  # Put info together to calculate GC
  sumv <- left + right + front + back + up + down
  local_gc <- sumv / sum_diff2 * (n - 1) / (2 * weights)
  global_gc <- sum(local_gc)

  return(global_gc)
}


# Calculate Geary's C in 2D
geary2D<-function(mat){

  # Set up
  n <- length(mat[is.na(mat)==F])
  xbar <- mean(mat, na.rm=T)
  diff <- mat - xbar
  diff2 <- diff**2
  sum_diff2 <- sum(diff2, na.rm = T)

  # Dimension of matrix
  nrow <- dim(diff)[1]
  ncol <- dim(diff)[2]

  # Rook shift
  up <- (mat - rbind(mat[-1,],rep(NA,ncol)) )^2
  down <- (mat - rbind(rep(NA,ncol),mat[-nrow,]) )^2
  left <- (mat - cbind(mat[,-1],rep(NA,nrow)) )^2
  right <- (mat - cbind(rep(NA,nrow),mat[,-ncol]) )^2

  # Calcultate weights
  weights <- sum(length(left[is.na(left) == F])) +
    sum(length(right[is.na(right) == F])) +
    sum(length(up[is.na(up) == F])) +
    sum(length(down[is.na(down) == F]))

  # Change all "NAs" to zero, so we can efficiently sum matrices
  left[is.na(left)] <- 0
  right[is.na(right)] <- 0
  up[is.na(up)] <- 0
  down[is.na(down)] <- 0

  # Put info together to calculate GC
  sumv <- left + right + up + down
  local_gc <- sumv / sum_diff2 * (n - 1) / (2 * weights)
  global_gc <- sum(local_gc)

  return(global_gc)
}






# Calculate fractal dimension in 2D
fd2D <- function(data) {
  fdx<-apply(data,1,fd1)
  fdy<-apply(data,2,fd1)
  fd<-c(fdx,fdy)
  fd_sub<-fd[fd<2]
  FD <- ifelse(length(fd_sub)<100,NA,1 + median(fd_sub,na.rm=T))
  return(FD)
}
fd1<-function(i){
  n <- length(i)
  g1 <- abs(i[2:n]-i[1:(n-1)])
  sumg1<-sum(g1,na.rm=T)/(2*length(g1[is.na(g1)==F]))
  g2 <- abs(i[3:n]-i[1:(n-2)])
  sumg2 <-  sum(g2,na.rm=T)/(2*length(g2[is.na(g2)==F]))
  slope <- sum(log(sumg2)-log(sumg1),na.rm=T)/log(2)
  fd <- 2-slope
  return(fd)
}