R/plot.BrcParcellation.R

In cdgreenidge/brcvis:

#' Plot method for objects of class `"BrcParcellation"`
#'
#' Plots an evenly-spaced series of slices through the brain. Each parcel is
#' colored differently. Note that empty slices are trimmed before plotting.
#'
#' @param x an object of class `"BrcParcellation"`
#' @param numSlices the number of slices to plot. Values greater than or equal
#'   to the maximum number of slices will plot all the slices.
#' @param view one of three strings: "sagittal", "coronal", or "axial".
#'   Specifies whether to plot sagittal, coronal, or axial slices.
#' @param colors A list of R colors similar to that generated by
#'   grDevices::rainbow().  Must contain a background color, plus exactly as
#'   many colors as there are unique parcels. If NULL, a default coloring is
#'   used.
#' @param ... not used
#' @export
plot.BrcParcellation <- function(x, numSlices, view="sagittal", colors=NULL,
                                 ...) {
  numParcels <- .numParcels(x)

  tryCatch({ brcbase::isValid(x) }, error=function(e) {
    stop(paste("Tried to plot invalid BrcParcellation object: ", e))
  })

  if ((numSlices %% 1 != 0) || (numSlices < 0)) {
    stop("numSlices argument must be a positive integer")
  }

  views <- list(sagittal=1, coronal=2, axial=3)
  if (!(view %in% names(views))) {
    stop(c("view argument must be one of 'sagittal', 'coronal', or 'axial'"))
  }

  if (is.null(colors)) {
    colors <- .defaultColors(numParcels)
  } else if (!all(.isColor(colors))) {
    stop("color argument contains invalid colors")
  } else if (.numParcels(x) != (length(colors) - 1)) {
    stop(paste("colors argument must contain 1 more color than the number ",
               "of parcels in the parcellation"))
  }


  dimension <- views[[view]]

  arr <- .parcellationToArray(x)
  arr <- .removeZeroSlices(arr, dimension)
  indices <- .makeIndexSequence(max=dim(arr)[dimension], length=numSlices)
  slices <- .extractSlices(arr, indices, dimension)
  invisible(.plotSlices(slices, numParcels(x), colors))
}

.isColor <- function(colors) {
  unname(sapply(colors, function(x) {
    tryCatch(is.matrix(grDevices::col2rgb(x)), error=function(e) FALSE)
  }))
}

.defaultColors <- function(numParcels) {
  c("#000000FF", grDevices::rainbow(numParcels))
}

.parcellationToArray <- function(parcellation) {
  data <- .factorToNumeric(parcellation$partition)
  array(data=data, dim=parcellation$dim3d)
}

.factorToNumeric <- function(xs) {
  as.numeric(levels(xs))[xs]
}

.removeZeroSlices <- function(arr, dim) {
  nonzero <- apply(arr, dim, function(xs) any(xs != 0))

  # R is really really bad at handling high-dimensional arrays. In particular,
  # there is no way to subset along a dimension if we only know the dimension
  # at runtime, because the syntax varies by dimension: arr[nonzero, , ] for
  # the first dimension, arr[ , nonzero, ] for the second dimension, etc.
  # So, we have to construct the function call by hand. How lovely!

  args <- rep(list(bquote()), times=length(dim(arr)))
  args[[dim]] <- quote(nonzero)
  call <- as.call(c(as.name("["), quote(arr), args))

  # No, I'm not bitter

  eval(call)
}

.makeIndexSequence <- function(max, length) {
  round(seq(1, max, length.out=length))
}

.extractSlices <- function(arr, indices, dim) {
  .splitAlongDim(arr, dim)[indices]
}

.numParcels <- function(parcellation) {
  l <- levels(parcellation$partition)
  length(l[l != 0])
}

.plotSlices <- function(slices, numParcels, colors) {
  layout <- .plotLayout(numSlices=length(slices))
  maxParcel <- max(unlist(slices))
  graphics::par(mfrow=c(layout$nrow, layout$ncol), mar=rep(0.2, 4), bg="black")
  for (i in 1:length(slices)) {
    graphics::image(slices[[i]],
                    asp=ncol(slices[[i]]) / nrow(slices[[i]]),
                    breaks=(0:(maxParcel + 1)) - 0.5,
                    bty="n",
                    col=colors,
                    xaxt="n",
                    yaxt="n")
  }
}

.plotLayout <- function(numSlices) {
	nrow = ceiling(sqrt(numSlices / 2))
	ncol = ceiling(numSlices / nrow)
	list(nrow=nrow, ncol=ncol)
}

# This function borrowed from
# http://stackoverflow.com/questions/20198751/three-dimensional-array-to-list
# Thanks, internet!
.splitAlongDim <- function(arr, dim) {
  stats::setNames(lapply(split(arr, arrayInd(seq_along(arr), dim(arr))[, dim]),
                         array, dim=dim(arr)[-dim], dimnames(arr)[-dim]),
                  dimnames(arr)[[dim]])
}