\documentclass[a4paper]{article}
\usepackage{filecontents}
\begin{filecontents}{qmri.bib}
@article{Weiskopf2013,
Author = {Weiskopf, N. and Suckling, J. and Williams, G. and Correia, M. M. and Inkster, B. and Tait, R. and Ooi, C. and Bullmore, E. T. and Lutti, A.},
Doi = {10.3389/fnins.2013.00095},
Journal = {Front. Neurosci.},
Language = {eng},
Medline-Pst = {epublish},
Pages = {95},
Pmid = {23772204},
Title = {Quantitative multi-parameter mapping of {R1, PD(*), MT, and R2(*) at 3T}: a multi-center validation.},
Url = {https://dx.doi.org/10.3389/fnins.2013.00095},
Volume = {7},
Year = {2013},
Bdsk-Url-1 = {https://dx.doi.org/10.3389/fnins.2013.00095}}
@article{Weiskopf2014,
Author = {N. Weiskopf and M. F. Callaghan and O. Josephs and A. Lutti and S. Mohammadi},
Journal = {Front. Neurosci.},
Number = {278},
Pages = {1--10},
Title = {Estimating the apparent transverse relaxation time {(R2*)} from images with different contrasts {(ESTATICS)} reduces motion artifacts},
Volume = {8},
Year = {2014}}
@article{Tabelow2019,
Author = {K. Tabelow and E. Balteau and J. Ashburner and M. F. Callaghan and B. Draganski and G. Helms and F. Kherif and T. Leutritz and A. Lutti and Ch. Phillips and E. Reimer and L. Ruthotto and M. Seif and N. Weiskopf and G. Ziegler and S. Mohammadi},
Doi = {10.1016/j.neuroimage.2019.01.029},
Journal = {NeuroImage},
Keywords = {Quantitative MRI, In vivo histology, Microstructure, Multi-parameter mapping, Relaxometry, SPM toolbox},
Pages = {191 - 210},
Title = {hMRI -- A toolbox for quantitative MRI in neuroscience and clinical research},
Volume = {194},
Year = {2019},
Bdsk-Url-1 = {https://www.sciencedirect.com/science/article/pii/S1053811919300291},
Bdsk-Url-2 = {https://doi.org/10.1016/j.neuroimage.2019.01.029}}
@conference{Tabelow2017,
Author = {Tabelow, K. and D'Alonzo, Ch. and Ruthotto, L. and Callaghan, M. F. and Weiskopf, N. and Polzehl, J. and Mohammadi, S.},
Booktitle = {Proceedings of ISMRM 2017},
Date-Added = {2019-01-09 10:05:03 +0100},
Date-Modified = {2019-01-09 10:05:03 +0100},
Title = {Removing the estimation bias due to the noise floor in multi-parameter maps},
Year = {2017}}
@article{PoTa2016,
Author = {J. Polzehl and K. Tabelow},
Doi = {10.1080/01621459.2016.1222284},
Journal = {J. Am. Stat. Assoc.},
Number = {516},
Pages = {1480-1490},
Title = {Low {SNR} in diffusion {MRI} models},
Volume = {111},
Year = {2016}}
@article{Tabelow2014a,
Author = {K. Tabelow and H. U. Voss and J. Polzehl},
Institution = {WIAS-Berlin},
Journal = {Med. Image Anal.},
Owner = {tabelow},
Pages = {76--86},
Timestamp = {2014.10.01},
Title = {Local estimation of the noise level in {MRI} using structural adaptation},
Volume = {20},
Year = {2015}}
@techreport{Tabelow2017a,
Author = {Mohammadi, S. and D'Alonzo, Ch. and Ruthotto, L. and Polzehl, J. and Ellerbrock, I. and Callaghan, M. F. and Weiskopf, N. and Tabelow, K.},
Doi = {10.20347/WIAS.PREPRINT.2432},
Institution = {WIAS},
Number = {2432},
Title = {Simultaneous adaptive smoothing of relaxometry and quantitative magnetization transfer mapping},
Type = {Preprint},
Year = {2017},
Bdsk-Url-1 = {https://doi.org/10.20347/WIAS.PREPRINT.2432}}
@techreport{PoPaTa18,
Author = {J. Polzehl and K. Papafitsoros and K. Tabelow},
Doi = {10.20347/WIAS.PREPRINT.2520},
Institution = {WIAS},
Number = {2520},
Title = {Patch-wise adaptive weights smoothing},
Type = {Preprint},
Url = {https://www.wias-berlin.de/preprint/2432/wias\_preprints_2432.pdf},
Year = {2018},
Bdsk-Url-1 = {https://www.wias-berlin.de/preprint/2432/wias_preprints_2432.pdf},
Bdsk-Url-2 = {https://doi.org/10.20347/WIAS.PREPRINT.2520},
Bdsk-Url-3 = {https://www.wias-berlin.de/preprint/2432/wias%5C_preprints_2432.pdf}}
@Book{MRBIbook,
title = {Magnetic Resonance Brain Imaging: Modeling and Data Analysis Using R},
publisher = {Springer},
year = {2019},
author = {J\"org Polzehl and Karsten Tabelow},
series = {Use R!},
doi = {10.1007/978-3-030-29184-6}
}
\end{filecontents}
\usepackage[style=authoryear,backend=bibtex,url=false]{biblatex} %backend tells biblatex what you will be using to process the bibliography file
\addbibresource{qmri}
\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}\index{Packages!#1}}
\let\proglang=\textsf
\let\code=\texttt
%\VignetteIndexEntry{An example session for analyzing quantitative MRI data in the Multi-Parameter Mapping framework}
\title{An example session for analyzing quantitative MRI data in the Multi-Parameter Mapping framework}
\author{J\"org Polzehl and Karsten Tabelow}
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\setkeys{Gin}{width=\textwidth}
This document illustrates the workflow of analyzing quantitative Magnetic Resonance Imaging (qMRI)
data in the framework of Multi-Parameter Mapping (MPM) experiments~\parencite{Weiskopf2013}.
The example uses artificial noisy MPM data based on the analysis results of a publicly
available MPM dataset, see \url{https://hmri.info}: The original data was modeled
using the ESTATICS model~\parencite{Weiskopf2014}. Its four adaptively smoothed parameter
maps were used to generate the artificial data on a selected (small) sub cube with Rician
noise by creating corresponding signal echos following the ESTATICS model equation.
The artificial data is supplied with this package together with the smoothed quantitative
maps (relaxation rates $R_1$ and $R_2^\star$, proton density $P\!D$, and magnetization
transfer $M\!T$) from the original data for ground truth comparison.
Specifically, the data consist of, in total, 22 image volumes in NIfTI format that
correspond to multiple echos of three different imaging modalities, i.e., 8 echos
of $T_1$ weighted images, 8 echos of proton density ($P\!D$) weighted images and
6 images from a dual excitation FLASH magnetization transfer ($M\!T$) weighted
sequence. Each echo in the image volumes is ``acquired'' at an echo time, repetition time
and flip angle, that equal those of the original experiment. The quantitative maps
suffer from a $B_1$ transmit (and receive) bias, that can be corrected for, if the smooth
$B_1$ correction field that affects the local flip angle is known~\parencite{Tabelow2019}.
An estimate for this field is supplied with the package, too.
For an more extended introduction we refer to \cite{MRBIbook} Chapter 6.
\section{Reading the MPM data}
<<0,echo=FALSE>>=
old <- options(digits=3, warn=0)
on.exit(options(old))
@
First, we specify the directory where the data are stored within the package
<<1>>=
dataDir <- system.file("extdata", package = "qMRI")
@
The filenames of the data correspond to the weighting of the imaging. We create
the paths of the data files including the $B_1$ field map a mask file (which
in this cases describes the whole data cube, only):
<<2>>=
t1Names <- paste0("t1w_", 1:8, ".nii.gz")
mtNames <- paste0("mtw_", 1:6, ".nii.gz")
pdNames <- paste0("pdw_", 1:8, ".nii.gz")
t1Files <- file.path(dataDir, t1Names)
mtFiles <- file.path(dataDir, mtNames)
pdFiles <- file.path(dataDir, pdNames)
B1File <- file.path(dataDir, "B1map.nii.gz")
maskFile <- file.path(dataDir, "mask.nii.gz")
@
The acquisition parameters (echo time ($T\!E$) in milliseconds, repetition time
($T\!R$) in milliseconds and flip angle ($F\!A$) in degree) for each volume are
replicated from the original MPM data:
<<3>>=
TE <- c(2.3, 4.6, 6.9, 9.2, 11.5, 13.8, 16.1, 18.4,
2.3, 4.6, 6.9, 9.2, 11.5, 13.8,
2.3, 4.6, 6.9, 9.2, 11.5, 13.8, 16.1, 18.4)
TR <- rep(25, 22)
FA <- c(rep(21, 8), rep(6, 6), rep(6, 8))
@
The package \pkg{pkg} provides a function \code{readMPMData} to read the data
into an \proglang{R} session
<<4>>=
library(qMRI)
mpm <- readMPMData(t1Files, pdFiles, mtFiles,
maskFile,
TR = TR, TE = TE, FA = FA,
verbose = FALSE)
@
which creates an object \code{mpm} of class \code{"MPMdata"}.
\section{Parameter estimation (ESTATICS)}
The ESTATICS model \parencite{Weiskopf2014} describes the data
from an MPM sequence by the relaxation rate $R_2^\star$ for the
signal decay with the echo time $T\!E$ and the three absolute
values for $T\!E = 0$. Parameters in this model can be estimated by
<<5>>=
setCores(2)
modelMPM <- estimateESTATICS(mpm,
method = "NLR",
verbose = FALSE)
@
<<7,echo=FALSE>>=
# evaluate this here to be able to reduce the mask afterwards
setCores(2, reprt=FALSE)
modelMPMQLsp1 <- smoothESTATICS(modelMPM,
mpmData = extract(mpm, "ddata"),
kstar = 16,
alpha = 0.004,
patchsize = 1,
verbose = FALSE)
@
<<7a, echo=FALSE>>=
# reduce mask to save time
mask <- extract(mpm,"mask")
mask[,c(1:10,12:21),] <- FALSE
mpm <- qMRI:::setMPMmask(mpm, mask)
@
using nonlinear least-squares regression. However, the low
low signal-to-noise ratio (SNR) in the data leads to a bias for the
parameter estimates, see~\cite{PoTa2016} and \cite{Tabelow2017}.
Therefore the parameter estimation can alternatively be performed by
<<6>>=
sigma <- array(50, mpm$sdim)
modelMPMQL <- estimateESTATICS(mpm,
method = "QL",
sigma = sigma,
L = 1,
verbose = FALSE)
@
using a quasi-likelihood formulation~\parencite{PoTa2016, Tabelow2017}.
It avoids, in case of low SNR, the bias caused by the skewness of the
Rician signal distribution. The application of this method
requires to specify the scale parameter $\sigma$ of the Rician distribution,
which in this case of artificial data is explicitly known.
Alternatively, a map of scale parameters
could be estimated using function \code{awslsigmc} from package \pkg{dti}
<<eval=FALSE>>=
ddata <- extract(mpm,"ddata")
mask <- extract(mpm,"mask")
if(require(dti)) sigma <- awslsigmc(ddata[1,,,],16,mask)$sigma
@
see \cite{Tabelow2014a} for details of the method.
\section{Structural adaptive smoothing}
Data from MPM sequences suffer from noise hindering the quality of
the estimated quantitative maps. \cite{Tabelow2017a} introduced a
structural adaptive smoothing method that can be used to reduce the
variability of the estimated parameter maps and, if \code{mpmData}
is specified, the observed image data:
% this code has been evaluated before
<<7b,eval=FALSE>>=
# use setCores(ncores) to enable openMP parallelization
setCores(2,reprt=FALSE)
modelMPMQLsp1 <- smoothESTATICS(modelMPMQL,
mpmData = extract(mpm, "ddata"),
kstar = 16,
alpha = 0.004,
patchsize = 1,
verbose = FALSE)
@
The package \pkg{qMRI} implements an extension of the method that
uses the comparison of local patches for adaptive smoothing~\parencite{PoPaTa18}.
The resulting ESTATICS unsmoothed and smoothed parameter maps for the
central coronal slice can be illustrated by
<<8, fig = TRUE, width = 12, height = 6.5>>=
library(adimpro)
rimage.options(zquantiles = c(.01, .99), ylab = "z")
oldpar <- par(mfrow = c(2, 4),
mar = c(3, 3, 3, 1), mgp = c(2, 1, 0))
on.exit(par(oldpar))
pnames <- c("T1", "MT", "PD", "R2star")
for (i in 1:4) {
modelCoeff <- extract(modelMPMQL,"modelCoeff")
rimage(modelCoeff[i, , 11, ])
title(pnames[i])
}
for (i in 1:4) {
modelCoeff <- extract(modelMPMQLsp1,"modelCoeff")
rimage(modelCoeff[i, , 11, ])
title(paste("smoothed", pnames[i]))
}
@
The preceding code also applies the adaptive smoothing not only to the
parameter maps of the ESTATICS model, but also to the data itself.
From this new parameter maps can be estimated by the quasi-likelihood
formulation and thus avoiding the bias due to the skewness in the signal
distribution:
<<9a,echo=FALSE>>=
mpmsp1 <- mpm
ddata <- extract(modelMPMQLsp1,"smoothedData")
dim(ddata) <- c(dim(ddata)[1],prod(dim(ddata)[-1]))
mpmsp1$ddata <- ddata[,mpm$mask]
@
<<9>>=
modelMPM2 <- estimateESTATICS(mpmsp1,
method = "NLR",
L = 1,
verbose = FALSE)
@
Note, that we again employ the same Rician scale parameter, i.e.,
of the unprocessed data, cf. the discussion in \cite{PoTa2016}.
\section{Calculating quantitative maps}
Finally, from the (un)smoothed parameter maps of the ESTATICS model, we
can compute the quantitative $R_1$, $R_2^\star$, $P\!D$ and $M\!T$ maps
<<10>>=
qMRIMaps <- calculateQI(modelMPM,
b1File = B1File,
TR2 = 3.4)
qMRIQLMaps <- calculateQI(modelMPMQL,
b1File = B1File,
TR2 = 3.4)
qMRIQLSmoothedp1Maps <- calculateQI(modelMPMQLsp1,
b1File = B1File,
TR2 = 3.4)
qMRIMaps2 <- calculateQI(modelMPM2,
b1File = B1File,
TR2 = 3.4)
@
The calculation requires the $B_1$-bias field, if available, and the $T\!R_2$
parameter of the MTw sequence \parencite{Tabelow2019}.
We show the central coronal slice of the estimated quantitative maps together with the
maps for the ground truth used to generate the data
<<11, fig = TRUE, width = 12, height = 13>>=
library(oro.nifti)
zlim <- matrix(c(0, 0, 0, 3000,
1.5, 35, 2, 10000),
4, 2)
R1 <- readNIfTI(file.path(dataDir, "R1map.nii.gz"))
R2star <- readNIfTI(file.path(dataDir, "R2starmap.nii.gz"))
MT <- readNIfTI(file.path(dataDir, "MTmap.nii.gz"))
PD <- readNIfTI(file.path(dataDir, "PDmap.nii.gz"))
rimage.options(ylab = "z")
par(mfrow = c(4, 4),
mar = c(3, 3, 3, 1), mgp = c(2, 1, 0))
nmaps <- c("R1", "R2star", "MT", "PD")
rimage(R1[, 11, ], zlim = zlim[1, ],
main = paste("true", nmaps[1]))
rimage(R2star[, 11, ], zlim = zlim[2, ],
main = paste("true", nmaps[2]))
rimage(MT[, 11, ], zlim = zlim[3, ],
main = paste("true", nmaps[3]),
col = colMT)
rimage(PD[, 11, ], zlim = zlim[4, ],
main = paste("true", nmaps[4]))
qmap1 <- extract(qMRIQLMaps, nmaps)
for (i in 1:4) rimage(qmap1[[i]][, 11, ], zlim = zlim[i, ],
main = paste("Estimated", nmaps[i]),
col = if(i==3) colMT else grey(0:225/255))
qmap2 <- extract(qMRIQLSmoothedp1Maps, nmaps)
for (i in 1:4) rimage(qmap2[[i]][, 11, ], zlim = zlim[i, ],
main = paste("Smoothed", nmaps[i]),
col = if(i==3) colMT else grey(0:225/255))
qmap3 <- extract(qMRIMaps2, nmaps)
for (i in 1:4) rimage(qmap3[[i]][, 11, ], zlim = zlim[i, ],
main = paste("Smoothed data", nmaps[i]),
col = if(i==3) colMT else grey(0:225/255))
@
The color scheme \code{colMT} is designed to emphasize the grey matter (magenta)
/ white matter (yellow) contrast in the MT image with a \code{zlim = c(0,2)}.
Using the quasi-likelihood estimation method is, for low SNR, supposed to
reduce the bias caused by the skewness of the Rician distribution,
cf.~\cite{PoTa2016} and \cite{Tabelow2017}.
<<12>>=
qmap0 <- extract(qMRIMaps,nmaps)
mask <- extract(mpm,"mask")
cat("\n",
"Bias of NLR estimates\n",
"R1", mean((qmap0$R1-R1)[mask]),
"R2star", mean((qmap0$R2star-R2star)[mask]),
"MT", mean((qmap0$MT-MT)[mask]),
"PD", mean((qmap0$PD-PD)[mask]), "\n",
"Bias of QL estimates\n",
"R1", mean((qmap1$R1-R1)[mask]),
"R2star", mean((qmap1$R2star-R2star)[mask]),
"MT", mean((qmap1$MT-MT)[mask]),
"PD", mean((qmap1$PD-PD)[mask]), "\n")
@
Let's see which estimate performs best with respect to the root mean squared error (RMSE):
<<13>>=
cat("\n",
"Root mean squared error of NLR estimate\n",
"R1", sqrt(mean((qmap0$R1-R1)[mask]^2)),
"R2star", sqrt(mean((qmap0$R2star-R2star)[mask]^2)),
"MT", sqrt(mean((qmap0$MT-MT)[mask]^2)),
"PD", sqrt(mean((qmap0$PD-PD)[mask]^2)), "\n",
"Root mean squared error of QL estimate\n",
"R1", sqrt(mean((qmap1$R1-R1)[mask]^2)),
"R2star", sqrt(mean((qmap1$R2star-R2star)[mask]^2)),
"MT", sqrt(mean((qmap1$MT-MT)[mask]^2)),
"PD", sqrt(mean((qmap1$PD-PD)[mask]^2)),"\n",
"Root mean squared error of smoothed QL estimate\n",
"R1", sqrt(mean((qmap2$R1-R1)[mask]^2)),
"R2star", sqrt(mean((qmap2$R2star-R2star)[mask]^2)),
"MT", sqrt(mean((qmap2$MT-MT)[mask]^2)),
"PD", sqrt(mean((qmap2$PD-PD)[mask]^2)),"\n",
"Root mean squared error of estimate from smoothed data \n",
"R1", sqrt(mean((qmap3$R1-R1)[mask]^2)),
"R2star", sqrt(mean((qmap3$R2star-R2star)[mask]^2)),
"MT", sqrt(mean((qmap3$MT-MT)[mask]^2)),
"PD", sqrt(mean((qmap3$PD-PD)[mask]^2)),"\n")
@
For interpretation we need to compare this to the mean parameter values
<<14>>=
cat("Mean R1", mean(R1[mask]), "Mean R2star",
mean(R2star[mask]), "Mean MT", mean(MT[mask]),
"Mean PD", mean(PD[mask]),"\n")
@
We see here, that using the quasi-likelihood estimation is supposed to
reduce the bias of the estimated ESTATICS parameters, but does so
at the cost of a slightly increased variance. Structural adaptive
smoothing leads to a considerable reduction in RMSE.
Modelling of the spatially smoothed data has the additional effect of
reduced bias originating from data variability.
\printbibliography
\end{document}