R/aslDenoiseR.R

In neuroconductor-devel/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' WIP: data-driven denoising for ASL MRI
#'
#' Denoises regression based reconstruction of CBF from arterial spin labeling
#'
#' @param boldmatrix input bold matrix
#' @param targety target to predict
#' @param covariates motion or other parameters / nuisance variables
#' @param selectionthresh e.g. 0.1 take 10 percent worst variables for noise
#' estimation
#' @param maxnoisepreds integer search range e.g 1:10
#' @param polydegree eg 4 for polynomial nuisance variables or 'loess'
#' @param crossvalidationgroups prior defined or integer valued
#' @param scalemat boolean
#' @param noisepoolfun function to help select noise pool e.g. max
#' @param usecompcor boolean
#' @param verbose boolean
#' @return matrix is output
#' @author Avants BB
#' @examples
#'
#' # asl<-antsImageRead( getANTsRData("pcasl") )
#' set.seed(1)
#' nvox <- 10*10*10*20
#' dims <- c(10,10,10,20)
#' asl <- makeImage( dims , rnorm( nvox )+500 )
#' aslmean <- getAverageOfTimeSeries( asl )
#' aslmask <- getMask( aslmean  )
#' aslmat<-timeseries2matrix( asl, aslmask )
#' for ( i in 1:10 ) aslmat[,i*2]<-aslmat[,i*2]*2
#' asl<-matrix2timeseries( asl, aslmask, aslmat )
#' tc<-as.factor(rep(c("C","T"),nrow(aslmat)/2))
#' dv<-computeDVARS(aslmat)
#' dnz<-aslDenoiseR( aslmat, tc, covariates=dv, selectionthresh=0.1,
#'   maxnoisepreds=c(1:2), polydegree=2, crossvalidationgroups=2 )
#' testthat::expect_equal(dnz$R2atBestN, 7, tolerance = 0.5 )
#'
#' \dontrun{
#' # a classic regression approach to estimating perfusion
#' # not recommended, but shows the basic idea.
#' # see ?quantifyCBF for a better approach
#' perfmodel<-lm( aslmat ~ tc + dnz$noiseu  )
#' perfimg<-antsImageClone(aslmask)
#' perfimg[ aslmask == 1 ]<-bigLMStats( perfmodel )$beta[1,]
#' m0<-getAverageOfTimeSeries(asl)
#' ctl<-c(1:(nrow(aslmat)/2))*2
#' m0[ aslmask==1 ]<-colMeans(aslmat[ctl,])
#' pcasl.parameters<-list( sequence="pcasl", m0=m0 )
#' cbf <- quantifyCBF( perfimg, aslmask, pcasl.parameters )
#'
#' # default mode network example
#'
#' if ( ! exists("bold") )
#'   {
#'   bold = antsImageRead( getANTsRData("rsbold") )
#'   meanbold = getAverageOfTimeSeries( bold )
#'   boldmask = getMask(  meanbold )
#'   # map to mni
#'   mni = antsImageRead( getANTsRData("mni") )
#'   mniaal = antsImageRead( getANTsRData("mnia") )
#'   mymap = antsRegistration( meanbold * boldmask, mni, typeofTransform='SyNBold',
#'    verbose=1 )
#'   aalimg = antsApplyTransforms( meanbold, mniaal, mymap$fwdtransforms,
#'     interpolator='NearestNeighbor')
#'   data("aal",package="ANTsR")
#'   timeselect<-10:dim(bold)[4]
#'   if ( ! exists("moco") )
#'     {
#'     moco = antsMotionCalculation( bold, boldmask )
#'     }
#'   sbold = smoothImage( moco$moco_img , 3.0 )
#'   antsImageWrite( boldmask, 'boldmask.nii.gz' )
#'   antsImageWrite( meanbold, 'boldmean.nii.gz' )
#'   antsImageWrite( aalimg, 'boldaal.nii.gz' )
#'   boldmask = boldmask * thresholdImage( aalimg, 1, Inf )
#'   }
#' postcing<-aal$label_num[ grep( "Cingulum_Post", aal$label_name ) ]
#' postCingMask = maskImage( boldmask, aalimg,
#'   level = as.numeric(postcing), binarize=T )
#' mpostCingMask= antsImageClone( postCingMask ) * 0
#' mpostCingMask[ postCingMask == 0 ] = 1
#' boldmat = timeseries2matrix( sbold, boldmask*mpostCingMask )
#' boldmat = boldmat[timeselect,]
#' boldmat = frequencyFilterfMRI( boldmat, tr=antsGetSpacing(bold)[4], opt='trig' )
#' dmnvec<-( timeseries2matrix(  sbold, postCingMask )[timeselect,] )
#' dmnvec = rowMeans(
#'   frequencyFilterfMRI( dmnvec, tr=antsGetSpacing(bold)[4], opt='trig' ) )
#' dmnmat = matrix( dmnvec, ncol=1)
#' mocpar = moco$moco_params[ timeselect , 3:14 ]
#' dnz<-aslDenoiseR( boldmat, dmnvec, covariates=mocpar, selectionthresh=0.2,
#'   maxnoisepreds=c(2:10),  polydegree='loess',
#'   crossvalidationgroups=8 )
#' boldmat<-timeseries2matrix(sbold, boldmask)
#' boldmat<-boldmat[timeselect,]
#' boldmat = frequencyFilterfMRI( boldmat, tr=antsGetSpacing(bold)[4], opt='trig' )
#' mdl<-bigLMStats( lm( boldmat ~ dmnvec + dnz$covariates + dnz$noiseu ), 0.001 )
#' betas<-mdl$beta.t[1,]
#' betaImg = makeImage( boldmask, betas )
#' antsImageWrite( betaImg, 'dmnBetas.nii.gz' )
#' # this should give default mode network around beta = 12
#' }
#'
#' @export aslDenoiseR
aslDenoiseR <- function(
  boldmatrix,
  targety,
  covariates = NA,
  selectionthresh = 0.1,
  maxnoisepreds = 2:12,
  polydegree = 'loess',
  crossvalidationgroups = 4,
  scalemat = F,
  noisepoolfun = max,
  usecompcor = F,
  verbose = F ) {
  nvox <- ncol(boldmatrix)
  groups <- crossvalidationgroups
  if (length(groups) == 1) {
    kfolds <- groups
    groups <- c()
    grouplength <- round(nrow(boldmatrix)/kfolds) - 1
    for (k in 1:kfolds) groups <- c(groups, rep(k, grouplength))
    groups <- c(rep(1, nrow(boldmatrix) - length(groups)), groups)
  }
  ##### identify low reproducibility voxels
  getnoisepool <- function(x, frac = selectionthresh) {
    xord <- sort(x)
    l <- round(length(x) * frac)
    val <- xord[l]
#    return( x < val )
    return(x < val & x < 0)
  }

  # overall description of the method 1. regressors include: design + trends +
  # noise-pool 2. find noise-pool by initial cross-validation without noise
  # regressors 3. cross-validate predictions using different numbers of noise
  # regressors 4. select best n for predictors from noise pool 5. return the noise
  # mask and the value for n make regressors
  timevals <- NA
  if (all(is.na(timevals)))
    timevals <- 1:nrow(boldmatrix)
  #
  if (is.numeric(polydegree)) {
    if (polydegree > 0) {
      p <- stats::poly(timevals, degree = polydegree)
      if (!all(is.na(covariates))) {
        covariates <- cbind(data.matrix(covariates), p)
      } else covariates <- p
    }
  } else if (polydegree == 'loess') {
    timevals <- 1:nrow(boldmatrix)
    mean.ts <- apply(boldmatrix, 1, mean)
    myloess <- loess(mean.ts ~ timevals)
    p <- myloess$fitted
    if (!all(is.na(covariates))) {
      covariates <- cbind(data.matrix(covariates), p)
    } else covariates <- p
  }
  rawboldmat <- data.matrix(boldmatrix)
  svdboldmat <- residuals( lm( rawboldmat ~ 1 + covariates ) )
  ################### now redo some work w/new hrf
  R2base <- crossvalidatedR2( svdboldmat, targety, groups,
    covariates = NA )
  R2base <- apply(R2base, FUN = noisepoolfun, MARGIN = 2)
  noisepool <- getnoisepool(R2base)
  if (all(noisepool == TRUE)) {
    print("all voxels meet your pvalthresh - try increasing the value")
    return(NA)
  }
  if (all(noisepool == FALSE)) {
    print("zero voxels meet your pvalthresh - try decreasing the value")
    return(NA)
  } else if ( verbose ) print(paste("Noise pool has nvoxels=", sum(noisepool)))
  if (scalemat)
    svdboldmat <- scale(svdboldmat)
  ##### should the denoising be done per group / run ?
  if (usecompcor)
    noiseu <- compcor(svdboldmat, max(maxnoisepreds))
  if (!usecompcor)
    noiseu <- svd(svdboldmat[, noisepool], nv = 0, nu = max(maxnoisepreds))$u
  R2summary <- rep(0, length(maxnoisepreds))
  ct <- 1
  for (i in maxnoisepreds) {
    svdboldmat <- residuals( lm( rawboldmat ~ 1 + covariates + noiseu[, 1:i]))
    R2 <- crossvalidatedR2( svdboldmat, targety, groups,
      covariates = NA )
    R2max <- apply(R2, FUN = max, MARGIN = 2)
    if (ct == 1)
      R2perNoiseLevel <- R2max else R2perNoiseLevel <- cbind(R2perNoiseLevel, R2max)
    R2pos <- R2max[R2max > 0]
    R2summary[ct] <- median(R2pos)
    if ( verbose ) print(paste("NoiseU:", i, "MeanRSqrd", R2summary[ct]))
    ct <- ct + 1
  }
  scl <- 0.95
  if (max(R2summary, na.rm = T) < 0)
    scl <- 1.05
  mxt <- scl * max(R2summary, na.rm = T)
  bestn <- maxnoisepreds[which(R2summary > mxt)[1]]
  if (ct == 2)
    R2final <- R2perNoiseLevel
  if (ct > 2)
    R2final <- R2perNoiseLevel[, bestn - min(maxnoisepreds) + 1]
  return(list(n = bestn, R2atBestN = R2summary[bestn], noisepool = noisepool, R2base = R2base,
    R2final = R2final, noiseu = noiseu[, 1:bestn],
    covariates = covariates))
}