#' Radiomic Calculation on CT Slices

#'

#' Calculate radiomic features on the each 2D slice of the whole 3D lung, right and left lungs separately

#'

#' @param img CT scan in ANTs image file format

#' @param mask Mask of CT scan in ANTs image file format

#' @param sides Choose to calculate radiomic features on the right and/or left lungs. Note: Right lung = 1, left lung = 2, non-lung = 0

#' @param plane One of: axial, coronal, sagittal

#' @param featuresFirst First level radiomic features to calculate

#' @param featuresSpatial Spatial radiomic features to calculate

#' @param tidy Logical. If true, outputs a tidy dataframe with results. If false, outputs nested loop.

#' @param reduce Logical. If true, reduces the dimensions of the scan based on extent of mask using reduce_scan.

#'

#' @return Radiomic values from every slice in both the right and left lungs

#' @export

#'

radiomics_slice <- function(img,

                           mask,

                           sides = c("right", "left"),

                           plane = 'axial',

                           featuresFirst = c('mean', 'sd', 'skew', 'kurtosis', 'min', 'q1', 'median', 'q3', 'max','energy', 'rms', 'uniformity', 'entropy'),

                           featuresSpatial = c('mi', 'gc', 'fd'),

                           tidy = TRUE,

                           reduce = TRUE){

  featuresMask <- lapply(sides, function(side){

    if(side == "right"){mv = 1}

    if(side == "left"){mv = 2}

    mask2 <- mask == mv

    # Reduce scan (optional)

    if(reduce == TRUE){

      red <- reduce_scan(img, mask2)

      img2 <- red$img

      mask2 <- red$mask

      rm(red)

      gc()

    }else{img2 <- img}

    # Put image in array format and remove non-mask values

    img2 <- as.array(img2)

    mask2 <- as.array(mask2)

    img2[mask2 != 1] <- NA

    # Which plane?

    if(plane == 'axial'){p = 3}

    if(plane == 'coronal'){p = 2}

    if(plane == 'sagittal'){p = 1}

    # Calculate features

    ndim <- dim(img2)[p]

    features <- apply(img2, p, function(x){

      npixels <- length(x[!is.na(x)])

      if(length(featuresFirst)>0){

        features1 <- radiomics_first(x, featuresFirst)

      }else(features1 <- NULL)

      if(length(featuresSpatial)>0){

        features2 <- radiomics_spatial(x, featuresSpatial)

      }else(features2 <- NULL)

      features <- c(features1, features2)

      features <- features[c(featuresFirst, featuresSpatial)]

      features <- c(npixels = npixels, features)

      return(features)

    })

    names(features) <- paste0('slic_num_', 1:ndim)

    return(features)

  })

  names(featuresMask) <- sides

  if(tidy == TRUE){

    # Make a nice little data frame to output

    test2 = NULL

    for(i in 1:length(featuresMask)){

      test <- do.call('rbind', featuresMask[[i]])

      test <- cbind.data.frame(lung = names(featuresMask)[i],

                               slice_number = names(featuresMask[[i]]),

                               test)

      test$slice_number <- gsub("slic_num_", "", test$slice_number)

      test <- as.data.frame(sapply(test, as.numeric))

      nslic <- dim(test)[1]

      test$slice_percent <- test$slice_number/nslic * 100

      test2 <- rbind(test2, test)

    }

    featuresMask <- test2

    rownames(featuresMask) <- c()

  }

  return(featuresMask)

}