R/aslAveraging.R

In ANTsX/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Average ASL tag-control pairs to estimate perfusion
#'
#' This function averages arterial spin labeling (ASL) functional MRI
#' tag-control image pairs to estimate perfusion.
#' @param asl input asl image
#' @param mask in which to calculate perfusion
#' @param tc vector indicating which images are tagged and which are controls.
#' Strongly recommended if using \code{regression} or \code{bayesian} methods.
#' @param nuisance nuisance covariates to include in regression
#' @param method method to use for computing average.  One of \code{sincSubtract},
#'  \code{simpleSubtract}, \code{cubicSubtract}, \code{surroundSubtract},
#' \code{regression}, or \code{bayesian}. See \code{Details}.
#' @param ... additional parameters to pass to ASL averaging functions.
#'   See \code{Details}.
#'
#' @details
#' Two major types of methods are available for ASL signal averaging:
#' \describe{
#'    \item{Subtraction}{: All the subtraction methods are based on subtracting
#'     the tag from control images and averaging the result.  \code{simple}
#'     subtracts adjacent tag and control images.  The other methods use
#'     interpolation to obtain a subtracted time-series.  Sinc subtraction may
#'     be marginally more accurate than cubic interpolation, but takes much
#'     longer.  Surround subtraction uses linear interpolation and is fast.}
#'    \item{Regression}{: Regression regresses the time-series on a vector
#'    of alternating tag-control dummy variables.  Regression can incorporate
#'    nuisance covariates.  Bayesian regression incorporates regularization in
#'    the regression to encourage all voxels of the same tissue type to have
#'    similar perfusion values.}
#'  }
#' For \code{bayesian}, two more arguments are required:
#' \describe{
#'   \item{segmentation}{: a segmentation image}
#'   \item{tissuelist}{: a list of tissue probability images}
#' }
#' These would be as output from \code{atropos}; see \code{Examples} for a
#' sample usage.
#'
#' @author Kandel BM, Avants BB
#' @examples
#' set.seed(1)
#' nvox <- 8 * 8 * 8 * 10
#' dims <- c(8, 8, 8, 10)
#' voxvals <- array(rnorm(nvox) + 500, dim = dims)
#' asl <- makeImage(dims, voxvals)
#' tc <- rep(c(-0.5, 0.5), dims[4] / 2)
#' avg <- aslAveraging(asl, tc = tc)
#' testthat::expect_equal(mean(avg), 0.00739, tolerance = .1)
#' slice <- extractSlice(asl, 4, 4)
#' mask <- getMask(slice)
#' seg <- atropos(d = 3, a = slice, x = mask, i = "kmeans[6]", m = "[0.0,1x1x1]")
#' perfSurr <- aslAveraging(asl, mask = NULL, method = "surroundSubtract")
#' bayesAvg <- aslAveraging(asl,
#'   tc = tc, method = "bayesian",
#'   segmentation = seg$segmentation, tissuelist = seg$probabilityimages
#' )
#'
#' @export
aslAveraging <- function(asl, mask = NULL, tc = NA, nuisance = NA, method = "regression", ...) {
  # define helper function
  bayesianPerfusion <- function(asl, xideal, nuisance, segmentation, tissuelist,
                                myPriorStrength = 30.0,
                                useDataDrivenMask = 3,
                                localweights = FALSE, priorBetas = NA) {
    segmentation <- check_ants(segmentation)
    mask <- thresholdImage(segmentation, 1, Inf)
    aslmat <- timeseries2matrix(asl, mask)
    perfdf <- data.frame(
      xideal = xideal,
      nuis = nuisance
    )
    perfdf <- as.matrix(perfdf[, !is.na(colMeans(perfdf))])
    perfmodel <- lm(aslmat ~ perfdf)
    getpriors <- function(img, segmentation) {
      n <- max(segmentation)
      p <- rep(0, n)
      segvec <- (segmentation[segmentation > 0])
      for (i in 1:n) {
        p[i] <- median(img[segvec == as.numeric(i)])
      }
      return(p)
    }
    if (all(is.na(priorBetas))) {
      blm <- bigLMStats(perfmodel, includeIntercept = T)
      bayespriormatfull <- blm$beta
    } else {
      bayespriormatfull <- priorBetas
    }
    n <- max(segmentation) * nrow(bayespriormatfull)
    bayespriormat <- matrix(rep(0, n), nrow = max(segmentation))
    for (i in 1:ncol(bayespriormat)) {
      bayespriormat[, i] <- getpriors(bayespriormatfull[i, ], segmentation)
    }
    #   set 4 to equal 2 - dgm = gm
    bayespriormat[4, ] <- bayespriormat[2, ]
    #   set csf to zero perfusion
    bayespriormat[1, 2] <- 0
    X <- model.matrix(perfmodel)
    localtissuemat <- imageListToMatrix(tissuelist, mask)
    priorwt <- diag(ncol(bayespriormat)) * myPriorStrength
    if (ncol(priorwt) > 2) {
      priorwt[3:ncol(priorwt), 3:ncol(priorwt)] <- 0
    }
    bayesianperfusionloc <- localtissuemat * 0
    bayesianperfusionlocp <- localtissuemat * 0
    for (i in 1:ncol(aslmat)) {
      # here is where we get really bayesian
      # average over all tissue models ...
      localtissuemat[, i] <- abs(localtissuemat[, i]) / sum(abs(localtissuemat[, i]))
      for (segval in 1:max(segmentation)) {
        tissueprior <- localtissuemat[segval, i]
        localprior <- bayespriormat[segval, ]
        blm <- bayesianlm(X, aslmat[, i], localprior, priorwt,
          includeIntercept = T
        )
        locbeta <- blm$beta[2]
        bayesianperfusionloc[segval, i] <- locbeta
        bayesianperfusionlocp[segval, i] <- locbeta * tissueprior
      }
    }
    bperfimg <- makeImage(mask, colSums(bayesianperfusionlocp))
    bperfimg
  }

  if ((method == "regression" | method == "bayesian") & all(is.na(tc))) {
    warning(paste(
      "Using regression to estimate perfusion, but not provided",
      "with \n tag-control labels.  Assuming first image is tag, \n",
      "with alternating tag-control pairs."
    ))
    labelfirst <- TRUE
    if (!labelfirst) {
      tc <- (rep(
        c(1, 0),
        dim(asl)[4]
      )[1:dim(asl)[4]] - 0.5) # control minus tag
    } else {
      tc <- (rep(
        c(0, 1),
        dim(asl)[4]
      )[1:dim(asl)[4]] - 0.5) # tag minus control
    }
  }

  if (length(grep("Subtract", method)) > 0) {
    avg <- ANTsRCore::timeSeriesSubtraction(asl, method)
  } else if (method == "regression") {
    labelfirst <- TRUE
    if (is.null(mask)) {
      myar <- apply(as.array(asl), c(1, 2, 3), mean)
      img <- makeImage(dim(myar), myar)
      antsSetSpacing(img, antsGetSpacing(asl)[1:3])
      antsSetOrigin(img, antsGetOrigin(asl)[1:3])
      antsSetDirection(img, antsGetDirection(asl)[1:3, 1:3])
      mask <- getMask(img)
    }
    ts <- timeseries2matrix(asl, mask)
    xideal <- tc
    cbfform <- formula(ts ~ xideal)
    if (!all(is.na(nuisance))) {
      cbfform <- formula(ts ~ xideal + nuisance)
    }
    mycbfmodel <- lm(cbfform) # standard regression
    cbfi <- antsImageClone(mask)
    betaideal <- ((mycbfmodel$coeff)[2, ])
    if (mean(betaideal) < 0) {
      betaideal <- (betaideal) * (-1)
    }
    cbfi[mask == 1] <- betaideal # standard results
    avg <- antsImageClone(cbfi)
  } else if (method == "bayesian") {
    avg <- bayesianPerfusion(asl, tc, nuisance, ...)
  } else {
    stop("Unrecognized method.")
  }
  if (mean(avg) < 0) {
    avg <- -avg
  }
  avg
}