R/radiomics_partition.R

In muschellij2/lungct: Processing of Lung Computed Tomography ('CT') Scans

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