sizing <- function(fn) paste0(round(file.info(fn)$size/1024^2, 2), "M")
statFile <- function(label, fn) paste0(label, "'", fn, "' ", sizing(fn), "\n")
statNifti <- function(label, im)
{
paste0(
label, "nifti[",
paste0(dim(im), collapse=" x "),
"] Pixel dimensions: ",
paste0(attr(im, "pixdim"), c(attr(im, "pixunits")[1], attr(im, "pixunits")), collapse=" x "),
'\n'
)
}
statMatrix <- function(label, mm)
{
paste0(
label, "matrix[",
paste0(dim(mm), collapse=" x "),
"]\n"
)
}
statVector <- function(label, mm)
{
paste0(
label, "vector[",
length(mm),
"]\n"
)
}
statInner <- function(label, obj)
{
if(is.null(obj)) return("")
if(all(is.na(obj))) return(paste0(label, "NA"))
if(class(obj)[1] == "character") return(statFile(label, obj))
if(class(obj)[1] == "niftiImage") return(statNifti(label, obj))
if(class(obj)[1] == "matrix") return(statMatrix(label, obj))
if(class(obj)[1] == "matrix") return(statVector(label, obj))
paste0(label, "Unhandled Class(",class(obj)[1],")\n")
}
#' @export
summary.statMap <- function(object, ...)
{
cat(paste0(
"\nFormula: ", paste0(as.character(object$formulas[[2]]), collapse=''), paste0(as.character(object$formulas[[1]]), collapse=''), "\n",
"\nContents:\n",
statInner(" Stat: ", object$stat),
statInner(" Coef: ", object$coef),
statInner(" Sqrt Sigma: ", object$sqrtSigma),
statInner(" Mask: ", object$mask),
statInner(" Template: ", object$template),
" Robust: ", object$robust, "\n"
))
}
#' @export
#' @method print statMap
print.statMap <- function(x, ...)
{
summary(x, ...)
}
redyellow = colorRampPalette(c('red', 'yellow'), space='Lab')
bluecyan = colorRampPalette(c('blue', 'cyan'), space='Lab')
#' Create images of a statMap
#'
#' modified from oro.nifti:::image.nifti
#'
#' @export
#' @param x the statMap object to display images of
#' @param thresh A threshold to apply to the image, defaults to 2.32
#' @param index Any selected image indexes to use for the z. defaults to NULL
#' @param col a vector of colors to use for scaled intensities defaults to a grey scale.
#' @param colpos a vector of colors to use for positive values.
#' @param colneg a vector of colors to use for negative values.
#' @param plane the plane to display, can be axial, coronal or sagittal.
#' @param xlab a title for the x axis.
#' @param ylab a title for the u axis.
#' @param axes display axes, defaults to false
#' @param oma A vector of the form c(bottom, left, top, right) giving the size of the outer margins in lines of text. Default to 0.
#' @param mar A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot. Defaults to 0.
#' @param bg background color, defaults to black.
#' @param ... additional arguments passed to par
#' @importFrom grDevices gray
#' @importFrom graphics par
# modified from oro.nifti:::image.nifti
#' @export
image.statMap = function (x, thresh=2.32, index = NULL, col = gray(0:64/64), colpos=redyellow(64), colneg=bluecyan(64),
plane = c("axial", "coronal", "sagittal"), xlab = "", ylab = "", axes = FALSE, oma = rep(0, 4), mar = rep(0, 4), bg = "black", ...)
{
object <- x
# mask can't be empty in typical statMap object unless it's manually constructed
if(is.null(object$mask)) object$mask=object$template
if(is.null(object$template)) object$template=object$mask
x = if(is.character(object$template)) readNifti(object$template) else object$template
pixdim = RNifti::pixdim(x)
mask = if(is.character(object$mask)) readNifti(object$mask) else object$mask
stat = if(is.character(object$stat)){
readNifti(object$stat)
} else if (length(dim(object$stat))!= length(dim(mask)))
stat.statMap(object)
else object$stat
switch(plane[1], axial = {
aspect <- pixdim[3]/pixdim[2]
}, coronal = {
x <- aperm(x, c(1, 3, 2))
mask <- aperm(mask, c(1, 3, 2))
stat <- aperm(stat, c(1, 3, 2))
aspect <- pixdim[4]/pixdim[2]
}, sagittal = {
x <- aperm(x, c(2, 3, 1))
mask <- aperm(mask, c(2, 3, 1))
stat <- aperm(stat, c(2, 3, 1))
aspect <- pixdim[4]/pixdim[3]
}, stop(paste("Orthogonal plane", plane[1], "is not valid.")))
# permuted image dimensions
# crop image and get image dimensions
xinds = apply(x!=0, 1, any)
yinds = apply(x!=0, 2, any)
zinds = apply(x!=0, 3, any)
x = x[xinds,,]
x = x[,yinds,]
x = x[,,zinds]
stat = stat[xinds,,]
stat = stat[,yinds,]
stat = stat[,,zinds]
statneg = stat
statneg[ statneg> -thresh] = 0
statneg = abs(statneg)
stat[ stat<thresh ] = 0
imgdim = dim(x)
zlim = range(x, na.rm=TRUE)
maxstat = max(c(stat[ stat>0], thresh), na.rm=TRUE)
maxstatneg = max(c(statneg[ statneg>0], thresh), na.rm=TRUE)
breaks <- c(zlim[1], seq(zlim[1], zlim[2], length = length(col) - 1), zlim[2])
breakspos <- c(thresh, seq(thresh, maxstat, length = length(colpos)-1), maxstat)
breaksneg <- c(thresh, seq(thresh, maxstatneg, length = length(colneg)-1), maxstatneg)
if(is.null(index)) index = 1:imgdim[3]
oldpar <- par(no.readonly = TRUE)
par(mfrow = ceiling(rep(sqrt(length(index)), 2)), oma = oma, mar = mar, bg = bg)
for (z in index) {
# background image
graphics::image(1:imgdim[1], 1:imgdim[2], x[, , z], col = col,
breaks = breaks, asp = aspect, axes = axes,
xlab = xlab, ylab = ylab, ...)
# overlay positive
graphics::image(1:imgdim[1], 1:imgdim[2], stat[, , z], col = colpos,
breaks = breakspos, asp = aspect, axes = axes, add=TRUE,
xlab = xlab, ylab = ylab, ...)
# overlay negative
graphics::image(1:imgdim[1], 1:imgdim[2], statneg[, , z], col = colneg,
breaks = breaksneg, asp = aspect, axes = axes, add=TRUE,
xlab = xlab, ylab = ylab, ...)
}
par(oldpar)
invisible()
}
#' Write the statMap objects out
#'
#' Given a statMap object and a directory write the objects as stat.nii.gz, coef.nii.gz and sqrtSigma.nii.gz
#' @param x the statMap object to write out
#' @param outdir the directory to write into
#' @return a list of what was written
#' @export
write.statMap <- function(x,outdir)
{
statimg = file.path(outdir, 'stat.nii.gz')
coefimg = file.path(outdir, 'coef.nii.gz')
resimg = file.path(outdir, 'sqrtSigma.nii.gz')
if(is.character(x$stat)){
file.copy(x$stat, statimg)
file.copy(x$sqrtSigma, resimg)
file.copy(x$coef, coefimg)
} else {
cat('Writing output images.\n')
dir.create(outdir, showWarnings=FALSE, recursive=TRUE)
writeNifti(stat.statMap(x), statimg)
writeNifti(coef.statMap(x), coefimg)
cat('Writing sqrtSigma 4d image.\n')
# reshape res matrix into 4d image
# memory intensive
# overwriting res
x$sqrtSigma = lapply(1:ncol(x$sqrtSigma), function(ind){ x$mask[ x$mask!=0] = x$sqrtSigma[,ind]; x$mask} )
# combine into 4d array
x$sqrtSigma = do.call(abind::abind, c(x$sqrtSigma, list(along=length(dim(x$mask))+1)))
writeNifti(updateNifti(x$sqrtSigma, x$mask), resimg)
res = resimg
}
return(list(stat=statimg, coef=coefimg, sqrtSigma=resimg))
}
#' Gets a 4D niftiImage of the coefficient image from a statMap object
#'
#' Returns a statistical niftiImage object from a statMap object
#' @param x the statMap object to extract a coefficient niftiImage from
#' @return a niftiImage object of the coefficient image
#' @export
stat.statMap = function(x){
if(is.character(x$stat)){
stat = readNifti(x$stat)
} else {
# output 4D coefficient image
stat = x$mask
stat[ stat!=0 ] = x$stat
}
return(stat)
}
#' Gets a 4D niftiImage of the statistical image from a statMap object
#'
#' Returns a 4D coefficient niftiImage object from a statMap object
#' @param object the statMap object to extract a coefficient niftiImage from
#' @param ... additional arguments (ignored)
#' @return a niftiImage object of the coefficient image
#' @export
coef.statMap = function(object, ...){
# output 4D coefficient image
coef = do.call(abind, c(lapply(1:nrow(object$coef), function(coefv){ object$mask[object$mask!=0] = coefv; object$mask}), list('along'=(length(dim(object$mask))+1)) ))
coef = updateNifti(coef, template=object$mask)
return(coef)
}
#' Gets a niftiImage of the variance image from a statMap object
#'
#' Will return a niftiImage of the variance image
#' @param x the statMap to extract the variance image from
#' @return a niftiImage object of the variance image
#' @export
var.statMap = function(x){
if(x$df>1)
stop('Only supported for df<1.')
# get niftiImage from mask
varimg = x$mask
varimg[ varimg!=0] = rowSums(x$sqrtSigma^2)/x$rdf
return(varimg)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.