R/bayesianCBF.R

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

#' Uses probabilistic segmentation to constrain pcasl-based cbf computation.
#'
#' Employs a robust regression approach to learn the relationship between a
#' sample image and a list of images that are mapped to the same space as the
#' sample image.
#'
#'
#' @param pcasl img antsImage for cbf
#' @param segmentation image, should cover the brain.
#' @param tissuelist a list containing antsImages eg list(prob1,...,probN)
#' @param myPriorStrength - e.g 30
#' @param useDataDrivenMask - morphology parameters e.g. 3
#' @param denoisingComponents - data-driven denoising parameters
#' @param robustnessvalue - value (e.g. 0.95) that throws away time points
#' @param localweights Use estimate of voxel-wise reliability to inform prior
#' weight?
#' @param priorBetas prior betas for each tissue and predictor
#' @return estimated cbf image
#' @author Brian Beaumont Avants and Benjamin M. Kandel
#' @keywords asl, bayesian blood cerebral flow,
#' @examples
#' \dontrun{
#' set.seed(123)
#' # see fMRIANTs github repository for data and I/O suggestions
#' }
#'
#' @export bayesianCBF
bayesianCBF <- function(pcasl, segmentation, tissuelist,
                        myPriorStrength = 30.0,
                        useDataDrivenMask = 3,
                        denoisingComponents = 1:8,
                        robustnessvalue = 0.95, # higher rejects more data. 0.9 or less - keep all
                        localweights = FALSE, priorBetas = NA) {
  compcorComponents <- 0
  motionAcc <- 2 # motion accuracy - 0 is for testing, 1 or 2 real studies
  if (all(dim(tissuelist[[1]]) == 1) | all(dim(seg) == 1) | all(dim(pcasl) == 1)) {
    stop(paste("Check your input data"))
  }
  avg <- getAverageOfTimeSeries(pcasl)
  avg <- n3BiasFieldCorrection(avg, 2)
  avg <- n3BiasFieldCorrection(avg, 2)
  mask <- antsImageClone(seg)
  mask[mask > 0] <- 1
  if (useDataDrivenMask > 0) {
    mask2 <- getMask(avg, mean(avg), Inf, useDataDrivenMask)
    # cleans up mask to agree with data-driven mask2
    mask[mask2 == 0] <- 0
    seg[mask2 == 0] <- 0
  }
  aslmat <- timeseries2matrix(pcasl, mask)
  perfpro <- aslPerfusion(pcasl,
    skip = 10,
    dorobust = robustnessvalue, useDenoiser = denoisingComponents,
    moreaccurate = motionAcc, verbose = 1, mask = mask, useBayesian = 0,
    ncompcor = compcorComponents
  )
  perfpro$m0 <- n3BiasFieldCorrection(perfpro$m0, 2)
  pcasl.parameters <- list(sequence = "pcasl", m0 = perfpro$m0)
  perfdf <- data.frame(
    xideal = perfpro$xideal,
    nuis = perfpro$nuisancevariables
  )
  perfdf <- perfdf[, !is.na(colMeans(perfdf))]
  perfmodel <- lm(aslmat ~ ., data = perfdf, weights = perfpro$regweights)
  getpriors <- function(img, seg) {
    n <- max(seg)
    p <- rep(0, n)
    #  Use the median to be conservative.
    segvec <- (seg[seg > 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 = TRUE)
    bayespriormatfull <- blm$beta
  }
  if (!all(is.na(priorBetas))) {
    bayespriormatfull <- priorBetas
  }
  n <- max(seg) * nrow(bayespriormatfull)
  bayespriormat <- matrix(rep(0, n), nrow = max(seg))
  for (i in 1:ncol(bayespriormat)) {
    bayespriormat[, i] <- getpriors(bayespriormatfull[i, ], seg)
  }
  # 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
  priorwt[3:ncol(priorwt), 3:ncol(priorwt)] <- 0
  ## alternative priors below
  # instead use bayespriormatfull - to est cov?
  priorwt2 <- solve(cov(bayespriormat) +
    diag(ncol(bayespriormat)) * .1) * myPriorStrength
  bayesianperfusionloc <- localtissuemat * 0
  bayesianperfusionlocp <- localtissuemat * 0
  if (localweights) {
    motion_params <- .motion_correction(pcasl, moreaccurate = 1)$moco_params[, 1:4]
    reliability <- aslDenoiseR(aslmat, perfpro$xideal, motion_params,
      usecompcor = T
    )$R2final
    reliability[reliability < 0.05] <- 0.05
    unreliability <- log(min(reliability) / reliability)
    if (max(unreliability) <= 0) {
      unreliability <- unreliability - min(unreliability)
    }
    unreliability <- unreliability / max(unreliability)
  }
  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(seg))
    {
      tissueprior <- localtissuemat[segval, i]
      localprior <- bayespriormat[segval, ]
      if (!localweights) {
        blm <- bayesianlm(X, aslmat[, i], localprior, priorwt,
          regweights = perfpro$regweights,
          includeIntercept = T
        )
      } else {
        blm <- bayesianlm(X, aslmat[, i], localprior, priorwt * unreliability[i],
          regweights = perfpro$regweights,
          includeIntercept = T
        )
      }
      locbeta <- blm$beta[2]
      bayesianperfusionloc[segval, i] <- locbeta
      bayesianperfusionlocp[segval, i] <- locbeta * tissueprior
    }
  }
  bperfimg <- makeImage(mask, colSums(bayesianperfusionlocp))
  bcbf <- quantifyCBF(bperfimg, mask, pcasl.parameters)
  bcbf <- bcbf$meancbf
  bayespriormatfull[2, ] <- colSums(bayesianperfusionlocp)
  return(list(
    bcbf = bcbf, nz = sum(perfpro$regweights < 0.001),
    posteriorBetaSummary = bayespriormatfull
  ))
}