#' 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)

}