#' Radiomic Calculation on Partitioned Lung

#'

#' Calculate radiomic features on the partitioned 3D lung

#'

#' @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 featuresFirst First level radiomic features to calculate

#' @param featuresSpatial Spatial radiomic features to calculate

#' @param partition Matrix of x, y, and z coordinates for each partition from partition_lung. If null, partition_lung is called.

#' @param kernel_size (If partition is null) Size of the kernel, in voxel units of width, depth, and height. Must be c(3,3,3) or greater. Default: c(30,30,30)

#' @param kernel_stride (If partition is null) Stride (or spacing) between kernels, in voxel units, for width, depth, and height. If kernel_stride = kernel_size, the partitions are non-overlapping. If stride = c(1,1,1), then each voxel is returned.

#' @param threshold Number of non-missing voxels needed to calculate radiomic features in each partition.

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

#'

#' @return Values from selected features for both left and right lungs

#' @importFrom ANTsR maskImage

#' @export

radiomics_partition <- function(img,

                                mask,

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

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

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

                                partition = NULL,

                                kernel_size = c(30, 30, 30),

                                kernel_stride = c(30, 30, 30),

                                threshold = 1000,

                                tidy = TRUE) {

  # Get partition, if necessary

  if(is.null(partition)){

    partition = partition_lung(img,

                               kernel_size = kernel_size,

                               kernel_stride = kernel_stride,

                               centroid = TRUE)

  }

  # Calculate radiomic features on partitions within each mask value

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

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

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

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

    img2 <- as.array(img)

    mask2 <- as.array(mask)

    mask2 <- mask2 == mv

    img2[mask2 != 1] <- NA

    # Calculate n each partition

    features <- lapply(1:dim(partition)[1], function(i){

      # Grab partition

      x <- img2[partition$x1[i]:partition$xend[i],

                partition$y1[i]:partition$yend[i],

                partition$z1[i]:partition$zend[i]]

      # Find dimension of partition and number of non-null pixels

      dim_p <- dim(x)

      if(is.null(dim_p)){dim_p <- c(0,0,0)}

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

      # Only calculate radiomic features if partition fits criteria

      if(dim_p[1] > 2 & dim_p[2] > 2 & dim_p[3] > 2 & npixels >= threshold){

        # Calculate features

        if(length(featuresFirst)>0){

          features1 <- radiomics_first(x, featuresFirst)

        }else(features1 <- NULL)

        if(length(featuresSpatial)>0){

          features2 <- radiomics_spatial(x, featuresSpatial)

        }else(features2 <- NULL)

        # Put features together and only keep specified features

        features <- c(features1, features2)

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

        features <- c(npixels = npixels, features)

      }else(features <- NULL)

      return(features)

    })

    # Name partitions and remove NULL partitions

    names(features) <- paste0('partition', 1:dim(partition)[1])

    features[sapply(features, is.null)] <- NULL

    return(features)

  })

  names(featuresMask) <- sides

  if(tidy == TRUE){

    # Make a nice little data frame to output

    test2 = NULL

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

      # Reduce list to data frame

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

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

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

                               test)

      # Reformatting

      test$partition <- gsub("partition", "", test$partition)

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

      # Add in partition centroids

      partition2 <- partition[partition$partition %in% test$partition, 8:10]

      test <- cbind(partition2, test)

      test2 <- rbind(test2, test)

    }

    featuresMask <- test2

    rm(test,test2,partition2)

    gc()

    rownames(featuresMask) <- c()

  }

  return(featuresMask)

}