#' 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 ) )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.