In muschellij2/oro.nifti: Rigorous - 'NIfTI' + 'ANALYZE' + 'AFNI' : Input / Output
library("bitops")
library("XML")
library("splines")
library("oro.nifti")
options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE)
# options(niftiAuditTrail = TRUE)
mniRL_jpeg = file.path(tempdir(), "mniRL.jpeg")
mniLR_jpeg = file.path(tempdir(), "mniLR.jpeg")
ffd_jpeg = file.path(tempdir(), "ffd.jpeg")
ffd_orthographic_jpeg = file.path(tempdir(), "ffd_orthographic.jpeg")
ffd_design_jpeg = file.path(tempdir(), "ffd_design.jpeg")
ffd_zstat1_jpeg = file.path(tempdir(), "ffd_zstat1.jpeg")
ffd_zstat2_jpeg = file.path(tempdir(), "ffd_zstat2.jpeg")
mniRL_orthographic_jpeg = file.path(tempdir(), "mniRL_orthographic.jpeg")
Introduction
Medical imaging is well established in both the clinical and research
areas with numerous equipment manufacturers supplying a wide variety
of modalities. The ANALYZE format was developed at the Mayo Clinic
(in the 1990s) to store multidimensional biomedical images. It is
fundamentally different from the DICOM standard since it groups all
images from a single acquisition (typically three- or
four-dimensional) into a pair of binary files, one containing header
information and one containing the image information. The DICOM
standard groups the header and image information, typically a single
two-dimensional image, into a single file. Hence, a single
acquisition will contain multiple DICOM files but only a pair of
ANALYZE files.
The NIfTI format was developed in the early 2000s by the DFWG (Data
Format Working Group) in an effort to improve upon the ANALYZE format.
The resulting NIfTI-1 format adheres to the basic header/image
combination from the ANALYZE format, but allows the pair of files to
be combined into a single file and re-defines the header fields. In
addition, NIfTI extensions allow one to store additional information
(e.g., key acquisition parameters, experimental design) inside a NIfTI
file.
\begin{table}[tbp]
\begin{center}
\begin{tabular}{p{0.475\textwidth}p{0.425\textwidth}}
\hline
\multicolumn{2}{c}{\pkg{oro.nifti}}\
\hline
\code{afni}, \code{anlz}, \code{nifti} & Class constructors for AFNI, ANALYZE and NIfTI objects.\
\code{as(, "nifti")} & Coerce object into class \code{nifti}.\
\code{audit.trail}, \code{aux.file}, \code{descrip} & Extract or replace slots in specific header fields.\
\code{fmri2oro}, \code{oro2fmri} & Convert between \code{fmridata} (\pkg{fmri}) and \code{nifti} objects.\
\code{hotmetal}, \code{tim.colors} & Useful color tables for visualization.\
\code{image}, \code{orthographic}, \code{overlay} & Two-dimensional visualization methods.\
\code{is.afni}, \code{is.anlz}, \code{is.nifti} & Logical checks.\
\code{readAFNI}, \code{readANALYZE}, \code{readNIfTI} & Data input.\
\code{writeAFNI}, \code{writeANALYZE}, \code{writeNIfTI} & Data output.\
\hline
\end{tabular}
\end{center}
\caption{List of functions available in \pkg{oro.nifti}.
Functionality around the AFNI data format was recently added to
the \pkg{oro.nifti} package. Please visit
\url{http://afni.nimh.nih.gov/afni} for more information about the
AFNI data format.}
\label{tab:functions}
\end{table}
The material presented here provides users with a method of
interacting with ANALYZE and NIfTI files in \proglang{R} \citep{R}.
Real-world data sets, that are publicly available, are used to
illustrate the basic functionality of \pkg{oro.nifti}
\citep{whi-sch-tho:JSS}. It should be noted that \pkg{oro.nifti}
focuses on functions for data input/output and visualization.
\proglang{S}4 classes \texttt{nifti} and \texttt{anlz} are provided
for further statistical analysis in \proglang{R} without losing
contextual information from the original ANALYZE or NIfTI files.
Images in the metadata-rich DICOM format may be converted to NIfTI
semi-automatically using \pkg{oro.dicom} by utilizing as much
information from the DICOM files as possible. Basic visualization
functions, similar to those commonly used in the medical imaging
community, are provided for \texttt{nifti} and \texttt{anlz} objects.
Additionally, the \pkg{oro.nifti} package allows one to track every
operation on a \texttt{nifti} object in an \proglang{XML}-based audit
trail.
The \pkg{oro.nifti} package should appeal not only to \proglang{R}
package developers, but also to scientists and researchers who want to
interrogate medical imaging data using the statistical capabilities of
\proglang{R} without writing and validating their own basic data
input/output functionality. Table~\ref{tab:functions} lists the key
functions for \pkg{oro.nifti} and groups them according to common
functionality. An example of using statistical methodology in
\proglang{R} for the analysis of functional magnetic resonance imaging
(fMRI) data is given in section~\ref{fmri_example}. Packages already
available on CRAN that utilize \pkg{oro.nifti} include:
\pkg{cudaBayesreg} \citep{daSilva:JSS}, \pkg{dcemriS4}
\citep{whi-sch:JSS}, \pkg{dpmixsim} \citep{dpmixsim}, and \pkg{RNiftyReg}
\citep{RNiftyReg}.
\section[oro.nifti: NIfTI-1 data input/output in R]{\pkg{oro.nifti}: NIfTI-1 data input/output in \proglang{R}}
Although the industry standard for medical imaging data is DICOM,
another format has come to be heavily used in the image analysis
community. The ANALYZE format was originally developed in conjunction
with an image processing system (of the same name) at the Mayo
Foundation. A common version of the format, although not the most
recent, is called ANALYZE 7.5. A copy of the file ANALYZE75.pdf has
been included in \pkg{oro.nifti} (accessed via
\code{system.file("doc/ANALYZE75.pdf", package="oro.dicom")}) since it
does not appear to be available from \url{www.mayo.edu} any longer.
An ANALYZE 7.5 format image is comprised of two files, the ".hdr"
and ".img" files, that contain information about the acquisition and
the acquisition itself, respectively. A more recent adaption of this
format is known as NIfTI-1 and is a product of the Data Format Working
Group (DFWG) from the Neuroimaging Informatics Technology Initiative
(NIfTI; \url{http://nifti.nimh.nih.gov}). The NIfTI-1 data format is
almost identical to the ANALYZE format, but offers a few improvements
\begin{itemize}
\item merging of the header and image information
into one file (.nii)
\item re-organization of the 348-byte fixed header into more
relevant categories
\item possibility of extending the header information.
\end{itemize}
There are several \proglang{R} packages that also offer input/output
functionality for the NIfTI and ANALYZE data formats in addition to
image analysis capabilities for specific MRI acquisition sequences;
e.g., \pkg{AnalyzeFMRI} \citep{AnalyzeFMRI}, \pkg{fmri}
\citep{pol-tab:fmri} and \pkg{tractor.base} \citep{tractor.base}. The
\pkg{Rniftilib} package provides access to NIfTI data via the
\proglang{nifticlib} library \citep{Rniftilib}.
\subsection{The NIfTI header}
\label{sec:nifti-header}
The NIfTI header inherits its structure (348 bytes in length) from the
ANALYZE data format. The last four bytes in the NIfTI header
correspond to the "magic" field and denote whether or not the header
and image are contained in a single file (\code{magic =
"n+1\textbackslash{}0"}) or two separate files (\code{magic =
"ni1\textbackslash{}0"}), the latter being identical to the
structure of the ANALYZE data format. The NIfTI data format added an
additional four bytes to allow for "extensions" to the header. By
default these four bytes are set to zero.
The first example of reading in, and displaying, medical imaging data
in NIfTI format \texttt{avg152T1_LR_nifti.nii.gz} was obtained from
the NIfTI website (\url{http://nifti.nimh.nih.gov/nifti-1/}).
Successful execution of the commands
fname <- system.file(file.path("nifti", "mniLR.nii.gz"), package="oro.nifti")
(mniLR <- readNIfTI(fname))
pixdim(mniLR)
descrip(mniLR)
aux.file(mniLR)
produces an \proglang{S}4 \code{"nifti"} object (or
\code{"niftiAuditTrail"} if the audit trail option is set). Some
accessor functions are also provided; e.g., \code{aux.file} and
\code{descrip}. The former is used to access the original name of the
file (if it has been provided) and the latter is the name of a valid
NIfTI header field used to hold a "description" (up to 80 characters
in length).
\subsection{The NIfTI image}
Image information begins at the byte position determined by the
\code{voxoffset} slot. In a single NIfTI file (\code{magic =
"n+1\textbackslash{}0"}), this is by default after the first 352
bytes. Header extensions extend the size of the header and come
before the image information leading to a consequent increase of
\code{voxoffset} for single NIfTI files. The split NIfTI (\code{magic
= "ni1\textbackslash{}0"}) and ANALYZE formats contain pairs of
files, where the header and image information are separated, and do
not have this problem. In this case \code{voxoffset} is set to 0.
The \code{image} function has been overloaded so that it behaves
differently when dealing with medical image objects (\code{nifti} and
\code{anlz}). The command
jpeg(filename=mniLR_jpeg, width=480, height=480, quality=95, bg="black")
image(mniLR)
dev.off()
image(mniLR)
produces a three-dimensional array of the MNI brain, with the default
NIfTI axes, and is displayed on a $10{\times}10$ grid of images
(Figure \ref{fig:mniLR+mniRL}a). The \code{image} function for
medical image \proglang{S}4 objects is an attempt to balance minimal
user input with enough flexibility to customize the display when
necessary. For example, single slices may be viewed by using the
option \code{plot.type="single"} in conjunction with the option
\code{z=} to specify the slice.
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.45\textwidth]{r mniLR_jpeg} &
\includegraphics[width=0.45\textwidth]{r mniRL_jpeg}\
\textbf{(a)} & \textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of MNI volume \code{mniLR_nifti}
stored in the \emph{neurological} convention (right-is-right), but
displayed in the \emph{radiological} convention (right-is-left).
\textbf{(b)} Axial slices of MNI volume \code{mniRL_nifti} stored
and displayed in the \emph{radiological} convention.}
\label{fig:mniLR+mniRL}
\end{figure}
The second example of reading in and displaying medical imaging data
in the NIfTI format \texttt{avg152T1_RL_nifti.nii.gz} was also
obtained from the NIfTI website
(\url{http://nifti.nimh.nih.gov/nifti-1/}). Successful execution of
the commands
fname <- system.file(file.path("nifti", "mniRL.nii.gz"), package="oro.nifti")
(mniRL <- readNIfTI(fname))
jpeg(filename=mniRL_jpeg, width=480, height=480, quality=95, bg="black")
image(mniRL)
dev.off()
image(mniRL)
produces a three-dimensional array of the MNI brain that is displayed
in a $10{\times}10$ grid of images (Figure \ref{fig:mniLR+mniRL}b).
The two sets of data in Figure \ref{fig:mniLR+mniRL} are stored in two
different orientations, commonly referred to as the
\emph{radiological} and \emph{neurological} conventions. The
neurological convention is where "right is right" and one is
essentially looking through the subject. The radiological convention
is where "right is left" and one is looking at the subject.
\begin{figure}[tbp]
\begin{center}
\includegraphics*[width=0.65\textwidth]{r mniRL_orthographic_jpeg}
\end{center}
\caption{Orthographic display of the MNI volume \code{mniRL_nifti}.
By default the mid-axial, mid-sagittal and mid-coronal planes are
chosen.}
\label{fig:mniRL-orthographic}
\end{figure}
An additional graphical display function has been added for
\code{nifti} and \code{anlz} objects that allows a so-called
orthographic visualization of the data.
jpeg(filename=mniRL_orthographic_jpeg, width=480, height=480, quality=95, bg="black")
orthographic(mniRL)
dev.off()
orthographic(mniRL)
As seen in Figure \ref{fig:mniRL-orthographic} the mid-axial,
mid-sagittal and mid-coronal planes are displayed by default. The
slices used may be set using \code{xyz = c(I,J,K)}, where $(I,J,K)$
are appropriate indices, and the crosshairs will provide a spatial
reference in each plane relative to the other two.
\subsection{A note on axes and orientation}
The NIfTI format contains an implicit generalized spatial
transformation from the data co-ordinate system $(i,j,k)$ into a
real-space "right-handed" co-ordinate system. In this real-space
system, the $(x,y,z)$ axes are \emph{usually} set such that $x$
increases from left to right, $y$ increases from posterior to anterior
and $z$ increases from inferior to superior.
At this point in time the \pkg{oro.nifti} package cannot apply an
arbitrary transform to the imaging data into $(x,y,z)$ space -- such a
transform may require non-integral indices and interpolation steps.
The package does accommodate straightforward transformations of
imaging data; e.g., setting the $i$-axis to increase from right to
left (the neurological convention). Future versions of
\pkg{oro.nifti} will attempt to address more complicated spatial
transformations and provide functionality to display the $(x,y,z)$
axes on orthographic plots.
\subsection[NIfTI and ANALYZE data in S4]{NIfTI and ANALYZE data in \proglang{S}4}
A major improvement in the \pkg{oro.nifti} package is the fact that
standard medical imaging formats are stored in unique classes under
the \proglang{S}4 system \citep{chambers:2008}. Essentially, NIfTI
and ANALYZE data are stored as multi-dimensional arrays with extra
slots created that capture the format-specific header information;
e.g., for a \code{nifti} object
slotNames(mniRL)
c(cal.min(mniRL), cal.max(mniRL))
range(mniRL)
mniRL@"datatype"
convert.datatype(mniRL@"datatype")
Note, an ANALYZE object has a slightly different set of slots. Slots
4--47 are taken verbatim from the definition of the NIfTI format and
are read directly from a file. The slot \code{.Data} is the
multidimensional array (since class \code{nifti} inherits from class
\code{array}) and the slots \code{trail}, \code{extensions} and
\code{reoriented} are used for internal bookkeeping. In the code
above we have accessed the min/max values of the imaging data using
the \code{cal.min} and \code{cal.max} accessor functions which matches
a direct interrogation of the \code{.Data} slot using the \code{range}
function. Looking at the \code{datatype} slot provides a numeric code
that may be converted into a value that indicates the type of byte
structure used (in this case an 8-bit or 1-byte unsigned integer).
As introduced in Section~\ref{sec:nifti-header} there are currently
only two accessor functions to slots in the NIfTI header
(\code{aux.file} and \code{descrip}) -- all other slots are either
ignored or used inside of functions that operate on ANALYZE/NIfTI
objects. The NIfTI class also has the ability to read and write
extensions that conform to the NIfTI data format. Customized printing
and validity-checking functions are available to the user and every
attempt has been made to ensure that the information from the
multi-dimensional array is in agreement with the header values.
The constructor function \code{nifti} produces valid NIfTI objects,
including a consistent header, from an arbitrary array.
n <- 100
(random.image <- nifti(array(runif(n*n), c(n,n,1))))
random.image@"dim_"
dim(random.image)
The function \code{writeNIfTI} outputs valid NIfTI class files, which
can be opened in other medical imaging software. Files can either be
stored as standard \code{.nii} files or compressed with gnuzip
(default).
writeNIfTI(random.image, "random")
list.files(pattern="random")
file.remove("random.nii.gz")
\subsection{The audit trail}
Following on from the \proglang{S}4 implementation of both the NIfTI
and ANALYZE data formats, the ability to extend the NIfTI data format
header is utilized in the \pkg{oro.nifti} package. Please use the
command
options(niftiAuditTrail=TRUE)
to turn on the "audit trail" option in \pkg{oro.nifti} and then
execute the function \code{enableAuditTrail()}. With the option
enabled extensions are properly handled when reading and writing NIfTI
data, users are allowed to add extensions to newly-created NIfTI
objects by casting them as \code{niftiExtension} objects and adding
\code{niftiExtensionSection} objects to the \code{extensions} slot,
and all operations that are performed on a NIfTI object will generate
what we call an \emph{audit trail} that consists of an
\proglang{XML}-based log \citep{XML}.
Figure \ref{fig:mniLR-audit-trail} displays output
from the accessor function \code{audit.trail(mniLR)}, the
\proglang{XML}-based audit trail that is stored as a NIfTI header
extension.
\begin{sidewaysfigure}
\centering
\begin{tabular}{p{22cm}}
audit.trail(mniLR)
\end{tabular}
\caption{\proglang{XML}-based audit trail obtained via
\code{audit.trail(mniLR)}. Note, this function will return
\code{NULL} if the \pkg{XML} package is not available.}
\label{fig:mniLR-audit-trail}
\end{sidewaysfigure}
Each log entry contains information not only about the function
applied to the NIfTI object, but also various system-level
information; e.g., version of \proglang{R}, user name, date, time,
etc. When writing NIfTI-class objects to disk, the
\proglang{XML}-based NIfTI extension is converted into plain text and
saved using \code{ecode=6} to denote plain ASCII text. The user may
control the tracking of data manipulation via the audit trail using the
global option \code{niftiAuditTrail}.
\subsection{Interactive visualization}
Basic visualization of \code{nifti} and \code{anlz} class images can
be achieved with any visualization for arrays in \proglang{R}. For
example, the \pkg{EBImage} package provides functions \code{display}
and \code{animate} for visualization \citep{EBImage}. Please note
that functions in \pkg{EBImage} expect grey-scale values in the range
$[0,1]$, hence the display of \code{nifti} data may be performed using
mniLR.range <- range(mniLR)
EBImage::display((mniLR - min(mniLR)) / diff(mniLR.range))
Interactive visualization of multi-dimensional arrays, stored in NIfTI
or ANALYZE format, is however best performed outside of \proglang{R}
at this point in time. Popular viewers, especially for neuroimaging
data, include
\begin{itemize}
\item FSLView (\url{http://www.fmrib.ox.ac.uk/fsl/fslview/}),
\item MRIcron (\url{http://cabiatl.com/mricron/}),
\item ITKSnap (\url{http://www.itksnap.org}), and
\item VolView (\url{http://www.kitware.com/products/volview.html}).
\end{itemize}
The \pkg{mritc} package provides basic interactive visualization of
ANALYZE/NIfTI data using a \proglang{Tcl}/\proglang{Tk} interface
\citep{mritc}.
\subsection{An example using functional MRI data}
\label{fmri_example}
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.6\textwidth]{r ffd_jpeg}\
\textbf{(a)}\
\includegraphics[width=0.6\textwidth]{r ffd_orthographic_jpeg}\
\textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of the functional MRI data set
\code{filtered_func_data} from the first acquisition.
\textbf{(b)} Orthographic display of the first volume from the
functional MRI data set \texttt{filtered_func_data}.}
\label{fig:ffd+orthographic}
\end{figure}
This is an example of reading in, and displaying, a four-dimensional
medical imaging data set in NIfTI format \texttt{filtered_func_data}
obtained from the \pkg{FSL} evaluation and example data suite
(\url{http://www.fmrib.ox.ac.uk/fsl/fsl/feeds.html}). Successful
execution of the commands
filtered.func.data <-
system.file(file.path("nifti", "filtered_func_data.nii.gz"),
package="oro.nifti")
(ffd <- readNIfTI(filtered.func.data))
jpeg(filename=ffd_jpeg, width=480, height=480, quality=95, bg="black")
image(ffd, zlim=range(ffd)*0.95)
dev.off()
image(ffd, zlim=range(ffd)*0.95)
produces a four-dimensional (4D) array of imaging data that may be
displayed in a $5{\times}5$ grid of images
(Figure \ref{fig:ffd+orthographic}a). The first three dimensions are
spatial locations of the voxel (volume element) and the fourth
dimension is time for this functional MRI (fMRI) acquisition. As seen
from the summary of object, there are 21 axial slices of fairly coarse
resolution ($4{\times}4{\times}6\;\text{mm}$) and reasonable temporal
resolution ($3\;\text{s}$). Figure \ref{fig:ffd+orthographic}b
depicts the orthographic display of the \texttt{filtered_func_data}
using the axial plane containing the left-and-right thalamus to
approximately center the crosshair vertically.
jpeg(filename=ffd_orthographic_jpeg, width=480, height=480, quality=95, bg="black")
orthographic(ffd, xyz=c(34,29,10), zlim=range(ffd)*0.9)
dev.off()
orthographic(ffd, xyz=c(34,29,10), zlim=range(ffd)*0.9)
\subsubsection{Statistical analysis}
The \proglang{R} programming environment provides a wide variety of
statistical methodology for the quantitative analysis of medical
imaging data. For example, functional MRI (fMRI) data are typically
analyzed by applying a multiple linear regression model, commonly
referred to in the literature as a general linear model (GLM), that
utilizes the stimulus experiment to construct the design matrix.
Estimation of the regression coefficients in the GLM produces a
statistical image; e.g., $Z$-statistics for a voxel-wise hypothesis
test on activation in fMRI experiments
\citep{fri-etal:spms,fri-etal:revisited}.
\begin{figure}[tbp]
\begin{center}
\includegraphics*[width=\textwidth]{r ffd_design_jpeg}
\end{center}
\caption{Visual (30 seconds on/off) and auditory (45 seconds on/off)
stimuli, convolved with a parametric haemodynamic response
function, used in the GLM-based fMRI analysis.}
\label{fig:ffd-design}
\end{figure}
The 4D volume of imaging data in \texttt{filtered_func_data} was
acquired in an experiment with a repetition time
$\text{TR}=3\;\text{s}$, using both visual and auditory stimuli. The
visual stimulus was applied using an on/off pattern for a duration of
60 seconds and the auditory stimulus was applied using an on/off
pattern for a duration of 90 seconds. A parametric haemodynamic
response function (HRF), with mean $\mu=6$ and standard deviation
$\sigma=3$, is utilized here which is similar to the default values in
\pkg{FSL} \citep{smi-etal:FSL}. We construct the experimental design and
HRF in seconds, perform the convolution and then downsample by a
factor of three in order to obtain columns of the design matrix that
match the acquisition of the MRI data.
visual <- rep(c(-0.5,0.5), each=30, times=9)
auditory <- rep(c(-0.5,0.5), each=45, times=6)
hrf <- c(dgamma(1:15, 4, scale=1.5))
hrf0 <- c(hrf, rep(0, length(visual)-length(hrf)))
visual.hrf <- convolve(hrf0, visual)
hrf0 <- c(hrf, rep(0, length(auditory)-length(hrf)))
auditory.hrf <- convolve(hrf0, auditory)
index <- seq(3, 540, by=3)
visual.hrf <- visual.hrf[index]
auditory.hrf <- auditory.hrf[index]
jpeg(ffd_design_jpeg, width=3*480, height=1.5*480, quality=95)
par(mfrow=c(1,2), mar=c(5,4,4,2) + 1, mex=0.85,
cex=1.5, cex.axis=1.5, cex.lab=1.5, cex.main=1.5)
plot(index, visual.hrf, type="l", lwd=2,
xlab="Acquisition Index", ylab="Visual Stimulus")
plot(index, auditory.hrf, type="l", lwd=2,
xlab="Acquisition Index", ylab="Auditory Stimulus")
dev.off()
Figure \ref{fig:ffd-design} depicts the visual and auditory stimuli,
convolved with the HRF, in the order of acquisition. The design
matrix is then used in a voxel-wise GLM, where the \code{lsfit}
function in \proglang{R} estimates the parameters in the linear
regression. At each voxel $t$-statistics and their associated
$p$-values are computed for the hypothesis test of no effect for each
individual stimulus, along with an $F$-statistic for the hypothesis
test of no effect of any stimuli using the \code{ls.print} function.
##reduced length due to R package storage limitations
visual.hrf<-visual.hrf[1:64]
auditory.hrf<-auditory.hrf[1:64]
## background threshold: 10% max intensity
voxel.lsfit <- function(x, thresh) { # general linear model
## check against background threshold
if (max(x) < thresh) {
return(rep(NA, 5))
}
## glm
output <- lsfit(cbind(visual.hrf, auditory.hrf), x)
## extract t-statistic, p-values
output.t <- ls.print(output, print.it=FALSE)$coef.table[[1]][2:3,3:4]
output.f <- ls.print(output, print.it=FALSE)$summary[3]
c(output.t, as.numeric(output.f))
}
## apply local glm to each voxel
ffd.glm <- apply(ffd, 1:3, voxel.lsfit, thresh=0.1 * max(ffd))
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.6\textwidth]{r ffd_zstat1_jpeg}\
\textbf{(a)}\
\includegraphics[width=0.6\textwidth]{r ffd_zstat2_jpeg}\
\textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of the functional MRI data with
the statistical image from the visual stimulus overlayed.
\textbf{(b)} Axial slices of the functional MRI data with the
statistical image from the auditory stimulus overlayed. Both sets
of test statistics were thresholded at $Z\geq{5}$ for all voxel.}
\label{fig:zstat1+zstat2}
\end{figure}
Given the multidimensional array of output from the GLM fitting
procedure, the $t$-statistics are separated and converted into
$Z$-statistics to follow the convention used in \pkg{FSL}. For the
purposes of this example we have not applied any multiple comparisons
correction procedure and, instead, use a relatively large threshold of
$Z>5$ for visualization.
dof <- ntim(ffd) - 1
Z.visual <- nifti(qnorm(pt(ffd.glm[1,,,], dof, log.p=TRUE), log.p=TRUE),
datatype=16)
Z.auditory <- nifti(qnorm(pt(ffd.glm[2,,,], dof, log.p=TRUE), log.p=TRUE),
datatype=16)
jpeg(ffd_zstat1_jpeg, width=480, height=480, quality=95, bg="black")
yrange <- c(5, max(Z.visual, na.rm=TRUE))
overlay(ffd, ifelse(Z.visual > 5, Z.visual, NA),
zlim.x=range(ffd)*0.95, zlim.y=yrange)
dev.off()
yrange <- c(5, max(Z.visual, na.rm=TRUE))
overlay(ffd, ifelse(Z.visual > 5, Z.visual, NA),
zlim.x=range(ffd)*0.95, zlim.y=yrange)
jpeg(ffd_zstat2_jpeg, width=480, height=480, quality=95, bg="black")
yrange <- c(5, max(Z.auditory, na.rm=TRUE))
overlay(ffd, ifelse(Z.auditory > 5, Z.auditory, NA),
zlim.x=range(ffd)*0.95, zlim.y=yrange)
dev.off()
yrange <- c(5, max(Z.auditory, na.rm=TRUE))
overlay(ffd, ifelse(Z.auditory > 5, Z.auditory, NA),
zlim.x=range(ffd)*0.95, zlim.y=yrange)
Statistical images in neuroimaging are commonly displayed as an
overlay on top of a reference image (one of the dynamic acquisitions)
in order to provide anatomical context. The \code{overlay} command in
\pkg{oro.nifti} allows one to display the statistical image of
voxel-wise activations overlayed on one of the original EPI (echo
planar imaging) volumes acquired in the fMRI experiment. The 3D array
of $Z$-statistics for the visual and auditory tasks are overlayed on
the original data for "anatomical" reference in
Figure \ref{fig:zstat1+zstat2}. The $Z$-statistics that exceed the
threshold appear to match know neuroanatomy, where the visual cortex
in the occipital lobe shows activation under the visual stimulus
(Figure \ref{fig:zstat1+zstat2}a) and the primary auditory cortex in
the temporal lobe shows activation under the auditory stimulus
(Figure \ref{fig:zstat1+zstat2}b).
\section{Conclusion}
Medical image analysis depends on the efficient manipulation and
conversion of DICOM data. The \pkg{oro.nifti} package has been
developed to provide the user with a set of functions that mask as
many of the background details as possible while still providing
flexible and robust performance.
The future of medical image analysis in \proglang{R} will benefit from
a unified view of the imaging data standards: DICOM, NIfTI, ANALYZE,
AFNI, MINC, etc. The existence of a single package for handling
imaging data formats would facilitate interoperability between the
ever increasing number of \proglang{R} packages devoted to medical
image analysis. We do not assume that the data structures in
\pkg{oro.nifti} are best-suited for this purpose and we welcome an
open discussion around how best to provide this standardization to the
end user.
\section*{Acknowledgments}
The authors would like to thank the National Biomedical Imaging
Archive (NBIA), the National Cancer Institute (NCI), the National
Institute of Health (NIH) and all institutions that have contributed
medical imaging data to the public domain. The authors would also
like to thank K. Tabelow for providing functionality around the AFNI
data format. VS is supported by the German Research Council (DFG
SCHM 2747/1-1).
muschellij2/oro.nifti documentation built on Nov. 28, 2023, 5:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library("bitops") library("XML") library("splines") library("oro.nifti") options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) # options(niftiAuditTrail = TRUE) mniRL_jpeg = file.path(tempdir(), "mniRL.jpeg") mniLR_jpeg = file.path(tempdir(), "mniLR.jpeg") ffd_jpeg = file.path(tempdir(), "ffd.jpeg") ffd_orthographic_jpeg = file.path(tempdir(), "ffd_orthographic.jpeg") ffd_design_jpeg = file.path(tempdir(), "ffd_design.jpeg") ffd_zstat1_jpeg = file.path(tempdir(), "ffd_zstat1.jpeg") ffd_zstat2_jpeg = file.path(tempdir(), "ffd_zstat2.jpeg") mniRL_orthographic_jpeg = file.path(tempdir(), "mniRL_orthographic.jpeg")
Introduction
Medical imaging is well established in both the clinical and research areas with numerous equipment manufacturers supplying a wide variety of modalities. The ANALYZE format was developed at the Mayo Clinic (in the 1990s) to store multidimensional biomedical images. It is fundamentally different from the DICOM standard since it groups all images from a single acquisition (typically three- or four-dimensional) into a pair of binary files, one containing header information and one containing the image information. The DICOM standard groups the header and image information, typically a single two-dimensional image, into a single file. Hence, a single acquisition will contain multiple DICOM files but only a pair of ANALYZE files.
The NIfTI format was developed in the early 2000s by the DFWG (Data Format Working Group) in an effort to improve upon the ANALYZE format. The resulting NIfTI-1 format adheres to the basic header/image combination from the ANALYZE format, but allows the pair of files to be combined into a single file and re-defines the header fields. In addition, NIfTI extensions allow one to store additional information (e.g., key acquisition parameters, experimental design) inside a NIfTI file.
\begin{table}[tbp]
\begin{center}
\begin{tabular}{p{0.475\textwidth}p{0.425\textwidth}}
\hline
\multicolumn{2}{c}{\pkg{oro.nifti}}\
\hline
\code{afni}, \code{anlz}, \code{nifti} & Class constructors for AFNI, ANALYZE and NIfTI objects.\
\code{as(
The material presented here provides users with a method of interacting with ANALYZE and NIfTI files in \proglang{R} \citep{R}. Real-world data sets, that are publicly available, are used to illustrate the basic functionality of \pkg{oro.nifti} \citep{whi-sch-tho:JSS}. It should be noted that \pkg{oro.nifti} focuses on functions for data input/output and visualization. \proglang{S}4 classes \texttt{nifti} and \texttt{anlz} are provided for further statistical analysis in \proglang{R} without losing contextual information from the original ANALYZE or NIfTI files. Images in the metadata-rich DICOM format may be converted to NIfTI semi-automatically using \pkg{oro.dicom} by utilizing as much information from the DICOM files as possible. Basic visualization functions, similar to those commonly used in the medical imaging community, are provided for \texttt{nifti} and \texttt{anlz} objects. Additionally, the \pkg{oro.nifti} package allows one to track every operation on a \texttt{nifti} object in an \proglang{XML}-based audit trail.
The \pkg{oro.nifti} package should appeal not only to \proglang{R} package developers, but also to scientists and researchers who want to interrogate medical imaging data using the statistical capabilities of \proglang{R} without writing and validating their own basic data input/output functionality. Table~\ref{tab:functions} lists the key functions for \pkg{oro.nifti} and groups them according to common functionality. An example of using statistical methodology in \proglang{R} for the analysis of functional magnetic resonance imaging (fMRI) data is given in section~\ref{fmri_example}. Packages already available on CRAN that utilize \pkg{oro.nifti} include: \pkg{cudaBayesreg} \citep{daSilva:JSS}, \pkg{dcemriS4} \citep{whi-sch:JSS}, \pkg{dpmixsim} \citep{dpmixsim}, and \pkg{RNiftyReg} \citep{RNiftyReg}.
\section[oro.nifti: NIfTI-1 data input/output in R]{\pkg{oro.nifti}: NIfTI-1 data input/output in \proglang{R}}
Although the industry standard for medical imaging data is DICOM, another format has come to be heavily used in the image analysis community. The ANALYZE format was originally developed in conjunction with an image processing system (of the same name) at the Mayo Foundation. A common version of the format, although not the most recent, is called ANALYZE 7.5. A copy of the file ANALYZE75.pdf has been included in \pkg{oro.nifti} (accessed via \code{system.file("doc/ANALYZE75.pdf", package="oro.dicom")}) since it does not appear to be available from \url{www.mayo.edu} any longer. An ANALYZE 7.5 format image is comprised of two files, the ".hdr" and ".img" files, that contain information about the acquisition and the acquisition itself, respectively. A more recent adaption of this format is known as NIfTI-1 and is a product of the Data Format Working Group (DFWG) from the Neuroimaging Informatics Technology Initiative (NIfTI; \url{http://nifti.nimh.nih.gov}). The NIfTI-1 data format is almost identical to the ANALYZE format, but offers a few improvements \begin{itemize} \item merging of the header and image information into one file (.nii) \item re-organization of the 348-byte fixed header into more relevant categories \item possibility of extending the header information. \end{itemize}
There are several \proglang{R} packages that also offer input/output functionality for the NIfTI and ANALYZE data formats in addition to image analysis capabilities for specific MRI acquisition sequences; e.g., \pkg{AnalyzeFMRI} \citep{AnalyzeFMRI}, \pkg{fmri} \citep{pol-tab:fmri} and \pkg{tractor.base} \citep{tractor.base}. The \pkg{Rniftilib} package provides access to NIfTI data via the \proglang{nifticlib} library \citep{Rniftilib}.
\subsection{The NIfTI header} \label{sec:nifti-header}
The NIfTI header inherits its structure (348 bytes in length) from the ANALYZE data format. The last four bytes in the NIfTI header correspond to the "magic" field and denote whether or not the header and image are contained in a single file (\code{magic = "n+1\textbackslash{}0"}) or two separate files (\code{magic = "ni1\textbackslash{}0"}), the latter being identical to the structure of the ANALYZE data format. The NIfTI data format added an additional four bytes to allow for "extensions" to the header. By default these four bytes are set to zero.
The first example of reading in, and displaying, medical imaging data in NIfTI format \texttt{avg152T1_LR_nifti.nii.gz} was obtained from the NIfTI website (\url{http://nifti.nimh.nih.gov/nifti-1/}). Successful execution of the commands
fname <- system.file(file.path("nifti", "mniLR.nii.gz"), package="oro.nifti") (mniLR <- readNIfTI(fname)) pixdim(mniLR) descrip(mniLR) aux.file(mniLR)
produces an \proglang{S}4 \code{"nifti"} object (or \code{"niftiAuditTrail"} if the audit trail option is set). Some accessor functions are also provided; e.g., \code{aux.file} and \code{descrip}. The former is used to access the original name of the file (if it has been provided) and the latter is the name of a valid NIfTI header field used to hold a "description" (up to 80 characters in length).
\subsection{The NIfTI image}
Image information begins at the byte position determined by the \code{voxoffset} slot. In a single NIfTI file (\code{magic = "n+1\textbackslash{}0"}), this is by default after the first 352 bytes. Header extensions extend the size of the header and come before the image information leading to a consequent increase of \code{voxoffset} for single NIfTI files. The split NIfTI (\code{magic = "ni1\textbackslash{}0"}) and ANALYZE formats contain pairs of files, where the header and image information are separated, and do not have this problem. In this case \code{voxoffset} is set to 0.
The \code{image} function has been overloaded so that it behaves differently when dealing with medical image objects (\code{nifti} and \code{anlz}). The command
jpeg(filename=mniLR_jpeg, width=480, height=480, quality=95, bg="black") image(mniLR) dev.off()
image(mniLR)
produces a three-dimensional array of the MNI brain, with the default NIfTI axes, and is displayed on a $10{\times}10$ grid of images (Figure \ref{fig:mniLR+mniRL}a). The \code{image} function for medical image \proglang{S}4 objects is an attempt to balance minimal user input with enough flexibility to customize the display when necessary. For example, single slices may be viewed by using the option \code{plot.type="single"} in conjunction with the option \code{z=} to specify the slice.
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.45\textwidth]{r mniLR_jpeg} &
\includegraphics[width=0.45\textwidth]{r mniRL_jpeg}\
\textbf{(a)} & \textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of MNI volume \code{mniLR_nifti}
stored in the \emph{neurological} convention (right-is-right), but
displayed in the \emph{radiological} convention (right-is-left).
\textbf{(b)} Axial slices of MNI volume \code{mniRL_nifti} stored
and displayed in the \emph{radiological} convention.}
\label{fig:mniLR+mniRL}
\end{figure}
The second example of reading in and displaying medical imaging data in the NIfTI format \texttt{avg152T1_RL_nifti.nii.gz} was also obtained from the NIfTI website (\url{http://nifti.nimh.nih.gov/nifti-1/}). Successful execution of the commands
fname <- system.file(file.path("nifti", "mniRL.nii.gz"), package="oro.nifti") (mniRL <- readNIfTI(fname))
jpeg(filename=mniRL_jpeg, width=480, height=480, quality=95, bg="black") image(mniRL) dev.off()
image(mniRL)
produces a three-dimensional array of the MNI brain that is displayed in a $10{\times}10$ grid of images (Figure \ref{fig:mniLR+mniRL}b). The two sets of data in Figure \ref{fig:mniLR+mniRL} are stored in two different orientations, commonly referred to as the \emph{radiological} and \emph{neurological} conventions. The neurological convention is where "right is right" and one is essentially looking through the subject. The radiological convention is where "right is left" and one is looking at the subject.
\begin{figure}[tbp]
\begin{center}
\includegraphics*[width=0.65\textwidth]{r mniRL_orthographic_jpeg}
\end{center}
\caption{Orthographic display of the MNI volume \code{mniRL_nifti}.
By default the mid-axial, mid-sagittal and mid-coronal planes are
chosen.}
\label{fig:mniRL-orthographic}
\end{figure}
An additional graphical display function has been added for \code{nifti} and \code{anlz} objects that allows a so-called orthographic visualization of the data.
jpeg(filename=mniRL_orthographic_jpeg, width=480, height=480, quality=95, bg="black") orthographic(mniRL) dev.off()
orthographic(mniRL)
As seen in Figure \ref{fig:mniRL-orthographic} the mid-axial, mid-sagittal and mid-coronal planes are displayed by default. The slices used may be set using \code{xyz = c(I,J,K)}, where $(I,J,K)$ are appropriate indices, and the crosshairs will provide a spatial reference in each plane relative to the other two.
\subsection{A note on axes and orientation}
The NIfTI format contains an implicit generalized spatial transformation from the data co-ordinate system $(i,j,k)$ into a real-space "right-handed" co-ordinate system. In this real-space system, the $(x,y,z)$ axes are \emph{usually} set such that $x$ increases from left to right, $y$ increases from posterior to anterior and $z$ increases from inferior to superior.
At this point in time the \pkg{oro.nifti} package cannot apply an arbitrary transform to the imaging data into $(x,y,z)$ space -- such a transform may require non-integral indices and interpolation steps. The package does accommodate straightforward transformations of imaging data; e.g., setting the $i$-axis to increase from right to left (the neurological convention). Future versions of \pkg{oro.nifti} will attempt to address more complicated spatial transformations and provide functionality to display the $(x,y,z)$ axes on orthographic plots.
\subsection[NIfTI and ANALYZE data in S4]{NIfTI and ANALYZE data in \proglang{S}4}
A major improvement in the \pkg{oro.nifti} package is the fact that standard medical imaging formats are stored in unique classes under the \proglang{S}4 system \citep{chambers:2008}. Essentially, NIfTI and ANALYZE data are stored as multi-dimensional arrays with extra slots created that capture the format-specific header information; e.g., for a \code{nifti} object
slotNames(mniRL) c(cal.min(mniRL), cal.max(mniRL)) range(mniRL) mniRL@"datatype" convert.datatype(mniRL@"datatype")
Note, an ANALYZE object has a slightly different set of slots. Slots 4--47 are taken verbatim from the definition of the NIfTI format and are read directly from a file. The slot \code{.Data} is the multidimensional array (since class \code{nifti} inherits from class \code{array}) and the slots \code{trail}, \code{extensions} and \code{reoriented} are used for internal bookkeeping. In the code above we have accessed the min/max values of the imaging data using the \code{cal.min} and \code{cal.max} accessor functions which matches a direct interrogation of the \code{.Data} slot using the \code{range} function. Looking at the \code{datatype} slot provides a numeric code that may be converted into a value that indicates the type of byte structure used (in this case an 8-bit or 1-byte unsigned integer).
As introduced in Section~\ref{sec:nifti-header} there are currently only two accessor functions to slots in the NIfTI header (\code{aux.file} and \code{descrip}) -- all other slots are either ignored or used inside of functions that operate on ANALYZE/NIfTI objects. The NIfTI class also has the ability to read and write extensions that conform to the NIfTI data format. Customized printing and validity-checking functions are available to the user and every attempt has been made to ensure that the information from the multi-dimensional array is in agreement with the header values.
The constructor function \code{nifti} produces valid NIfTI objects, including a consistent header, from an arbitrary array.
n <- 100 (random.image <- nifti(array(runif(n*n), c(n,n,1)))) random.image@"dim_" dim(random.image)
The function \code{writeNIfTI} outputs valid NIfTI class files, which can be opened in other medical imaging software. Files can either be stored as standard \code{.nii} files or compressed with gnuzip (default).
writeNIfTI(random.image, "random") list.files(pattern="random")
file.remove("random.nii.gz")
\subsection{The audit trail}
Following on from the \proglang{S}4 implementation of both the NIfTI and ANALYZE data formats, the ability to extend the NIfTI data format header is utilized in the \pkg{oro.nifti} package. Please use the command
options(niftiAuditTrail=TRUE)
to turn on the "audit trail" option in \pkg{oro.nifti} and then execute the function \code{enableAuditTrail()}. With the option enabled extensions are properly handled when reading and writing NIfTI data, users are allowed to add extensions to newly-created NIfTI objects by casting them as \code{niftiExtension} objects and adding \code{niftiExtensionSection} objects to the \code{extensions} slot, and all operations that are performed on a NIfTI object will generate what we call an \emph{audit trail} that consists of an \proglang{XML}-based log \citep{XML}.
Figure \ref{fig:mniLR-audit-trail} displays output from the accessor function \code{audit.trail(mniLR)}, the \proglang{XML}-based audit trail that is stored as a NIfTI header extension.
\begin{sidewaysfigure} \centering \begin{tabular}{p{22cm}}
audit.trail(mniLR)
\end{tabular} \caption{\proglang{XML}-based audit trail obtained via \code{audit.trail(mniLR)}. Note, this function will return \code{NULL} if the \pkg{XML} package is not available.} \label{fig:mniLR-audit-trail} \end{sidewaysfigure}
Each log entry contains information not only about the function applied to the NIfTI object, but also various system-level information; e.g., version of \proglang{R}, user name, date, time, etc. When writing NIfTI-class objects to disk, the \proglang{XML}-based NIfTI extension is converted into plain text and saved using \code{ecode=6} to denote plain ASCII text. The user may control the tracking of data manipulation via the audit trail using the global option \code{niftiAuditTrail}.
\subsection{Interactive visualization}
Basic visualization of \code{nifti} and \code{anlz} class images can be achieved with any visualization for arrays in \proglang{R}. For example, the \pkg{EBImage} package provides functions \code{display} and \code{animate} for visualization \citep{EBImage}. Please note that functions in \pkg{EBImage} expect grey-scale values in the range $[0,1]$, hence the display of \code{nifti} data may be performed using
mniLR.range <- range(mniLR) EBImage::display((mniLR - min(mniLR)) / diff(mniLR.range))
Interactive visualization of multi-dimensional arrays, stored in NIfTI or ANALYZE format, is however best performed outside of \proglang{R} at this point in time. Popular viewers, especially for neuroimaging data, include \begin{itemize} \item FSLView (\url{http://www.fmrib.ox.ac.uk/fsl/fslview/}), \item MRIcron (\url{http://cabiatl.com/mricron/}), \item ITKSnap (\url{http://www.itksnap.org}), and \item VolView (\url{http://www.kitware.com/products/volview.html}). \end{itemize} The \pkg{mritc} package provides basic interactive visualization of ANALYZE/NIfTI data using a \proglang{Tcl}/\proglang{Tk} interface \citep{mritc}.
\subsection{An example using functional MRI data} \label{fmri_example}
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.6\textwidth]{r ffd_jpeg}\
\textbf{(a)}\
\includegraphics[width=0.6\textwidth]{r ffd_orthographic_jpeg}\
\textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of the functional MRI data set
\code{filtered_func_data} from the first acquisition.
\textbf{(b)} Orthographic display of the first volume from the
functional MRI data set \texttt{filtered_func_data}.}
\label{fig:ffd+orthographic}
\end{figure}
This is an example of reading in, and displaying, a four-dimensional medical imaging data set in NIfTI format \texttt{filtered_func_data} obtained from the \pkg{FSL} evaluation and example data suite (\url{http://www.fmrib.ox.ac.uk/fsl/fsl/feeds.html}). Successful execution of the commands
filtered.func.data <- system.file(file.path("nifti", "filtered_func_data.nii.gz"), package="oro.nifti") (ffd <- readNIfTI(filtered.func.data))
jpeg(filename=ffd_jpeg, width=480, height=480, quality=95, bg="black") image(ffd, zlim=range(ffd)*0.95) dev.off()
image(ffd, zlim=range(ffd)*0.95)
produces a four-dimensional (4D) array of imaging data that may be displayed in a $5{\times}5$ grid of images (Figure \ref{fig:ffd+orthographic}a). The first three dimensions are spatial locations of the voxel (volume element) and the fourth dimension is time for this functional MRI (fMRI) acquisition. As seen from the summary of object, there are 21 axial slices of fairly coarse resolution ($4{\times}4{\times}6\;\text{mm}$) and reasonable temporal resolution ($3\;\text{s}$). Figure \ref{fig:ffd+orthographic}b depicts the orthographic display of the \texttt{filtered_func_data} using the axial plane containing the left-and-right thalamus to approximately center the crosshair vertically.
jpeg(filename=ffd_orthographic_jpeg, width=480, height=480, quality=95, bg="black") orthographic(ffd, xyz=c(34,29,10), zlim=range(ffd)*0.9) dev.off()
orthographic(ffd, xyz=c(34,29,10), zlim=range(ffd)*0.9)
\subsubsection{Statistical analysis}
The \proglang{R} programming environment provides a wide variety of statistical methodology for the quantitative analysis of medical imaging data. For example, functional MRI (fMRI) data are typically analyzed by applying a multiple linear regression model, commonly referred to in the literature as a general linear model (GLM), that utilizes the stimulus experiment to construct the design matrix. Estimation of the regression coefficients in the GLM produces a statistical image; e.g., $Z$-statistics for a voxel-wise hypothesis test on activation in fMRI experiments \citep{fri-etal:spms,fri-etal:revisited}.
\begin{figure}[tbp]
\begin{center}
\includegraphics*[width=\textwidth]{r ffd_design_jpeg}
\end{center}
\caption{Visual (30 seconds on/off) and auditory (45 seconds on/off)
stimuli, convolved with a parametric haemodynamic response
function, used in the GLM-based fMRI analysis.}
\label{fig:ffd-design}
\end{figure}
The 4D volume of imaging data in \texttt{filtered_func_data} was acquired in an experiment with a repetition time $\text{TR}=3\;\text{s}$, using both visual and auditory stimuli. The visual stimulus was applied using an on/off pattern for a duration of 60 seconds and the auditory stimulus was applied using an on/off pattern for a duration of 90 seconds. A parametric haemodynamic response function (HRF), with mean $\mu=6$ and standard deviation $\sigma=3$, is utilized here which is similar to the default values in \pkg{FSL} \citep{smi-etal:FSL}. We construct the experimental design and HRF in seconds, perform the convolution and then downsample by a factor of three in order to obtain columns of the design matrix that match the acquisition of the MRI data.
visual <- rep(c(-0.5,0.5), each=30, times=9) auditory <- rep(c(-0.5,0.5), each=45, times=6) hrf <- c(dgamma(1:15, 4, scale=1.5)) hrf0 <- c(hrf, rep(0, length(visual)-length(hrf))) visual.hrf <- convolve(hrf0, visual) hrf0 <- c(hrf, rep(0, length(auditory)-length(hrf))) auditory.hrf <- convolve(hrf0, auditory) index <- seq(3, 540, by=3) visual.hrf <- visual.hrf[index] auditory.hrf <- auditory.hrf[index]
jpeg(ffd_design_jpeg, width=3*480, height=1.5*480, quality=95) par(mfrow=c(1,2), mar=c(5,4,4,2) + 1, mex=0.85, cex=1.5, cex.axis=1.5, cex.lab=1.5, cex.main=1.5) plot(index, visual.hrf, type="l", lwd=2, xlab="Acquisition Index", ylab="Visual Stimulus") plot(index, auditory.hrf, type="l", lwd=2, xlab="Acquisition Index", ylab="Auditory Stimulus") dev.off()
Figure \ref{fig:ffd-design} depicts the visual and auditory stimuli, convolved with the HRF, in the order of acquisition. The design matrix is then used in a voxel-wise GLM, where the \code{lsfit} function in \proglang{R} estimates the parameters in the linear regression. At each voxel $t$-statistics and their associated $p$-values are computed for the hypothesis test of no effect for each individual stimulus, along with an $F$-statistic for the hypothesis test of no effect of any stimuli using the \code{ls.print} function.
##reduced length due to R package storage limitations visual.hrf<-visual.hrf[1:64] auditory.hrf<-auditory.hrf[1:64] ## background threshold: 10% max intensity voxel.lsfit <- function(x, thresh) { # general linear model ## check against background threshold if (max(x) < thresh) { return(rep(NA, 5)) } ## glm output <- lsfit(cbind(visual.hrf, auditory.hrf), x) ## extract t-statistic, p-values output.t <- ls.print(output, print.it=FALSE)$coef.table[[1]][2:3,3:4] output.f <- ls.print(output, print.it=FALSE)$summary[3] c(output.t, as.numeric(output.f)) } ## apply local glm to each voxel ffd.glm <- apply(ffd, 1:3, voxel.lsfit, thresh=0.1 * max(ffd))
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.6\textwidth]{r ffd_zstat1_jpeg}\
\textbf{(a)}\
\includegraphics[width=0.6\textwidth]{r ffd_zstat2_jpeg}\
\textbf{(b)}
\end{tabular}
\end{center}
\caption{\textbf{(a)} Axial slices of the functional MRI data with
the statistical image from the visual stimulus overlayed.
\textbf{(b)} Axial slices of the functional MRI data with the
statistical image from the auditory stimulus overlayed. Both sets
of test statistics were thresholded at $Z\geq{5}$ for all voxel.}
\label{fig:zstat1+zstat2}
\end{figure}
Given the multidimensional array of output from the GLM fitting procedure, the $t$-statistics are separated and converted into $Z$-statistics to follow the convention used in \pkg{FSL}. For the purposes of this example we have not applied any multiple comparisons correction procedure and, instead, use a relatively large threshold of $Z>5$ for visualization.
dof <- ntim(ffd) - 1 Z.visual <- nifti(qnorm(pt(ffd.glm[1,,,], dof, log.p=TRUE), log.p=TRUE), datatype=16) Z.auditory <- nifti(qnorm(pt(ffd.glm[2,,,], dof, log.p=TRUE), log.p=TRUE), datatype=16)
jpeg(ffd_zstat1_jpeg, width=480, height=480, quality=95, bg="black") yrange <- c(5, max(Z.visual, na.rm=TRUE)) overlay(ffd, ifelse(Z.visual > 5, Z.visual, NA), zlim.x=range(ffd)*0.95, zlim.y=yrange) dev.off()
yrange <- c(5, max(Z.visual, na.rm=TRUE)) overlay(ffd, ifelse(Z.visual > 5, Z.visual, NA), zlim.x=range(ffd)*0.95, zlim.y=yrange)
jpeg(ffd_zstat2_jpeg, width=480, height=480, quality=95, bg="black") yrange <- c(5, max(Z.auditory, na.rm=TRUE)) overlay(ffd, ifelse(Z.auditory > 5, Z.auditory, NA), zlim.x=range(ffd)*0.95, zlim.y=yrange) dev.off()
yrange <- c(5, max(Z.auditory, na.rm=TRUE)) overlay(ffd, ifelse(Z.auditory > 5, Z.auditory, NA), zlim.x=range(ffd)*0.95, zlim.y=yrange)
Statistical images in neuroimaging are commonly displayed as an overlay on top of a reference image (one of the dynamic acquisitions) in order to provide anatomical context. The \code{overlay} command in \pkg{oro.nifti} allows one to display the statistical image of voxel-wise activations overlayed on one of the original EPI (echo planar imaging) volumes acquired in the fMRI experiment. The 3D array of $Z$-statistics for the visual and auditory tasks are overlayed on the original data for "anatomical" reference in Figure \ref{fig:zstat1+zstat2}. The $Z$-statistics that exceed the threshold appear to match know neuroanatomy, where the visual cortex in the occipital lobe shows activation under the visual stimulus (Figure \ref{fig:zstat1+zstat2}a) and the primary auditory cortex in the temporal lobe shows activation under the auditory stimulus (Figure \ref{fig:zstat1+zstat2}b).
\section{Conclusion}
Medical image analysis depends on the efficient manipulation and conversion of DICOM data. The \pkg{oro.nifti} package has been developed to provide the user with a set of functions that mask as many of the background details as possible while still providing flexible and robust performance.
The future of medical image analysis in \proglang{R} will benefit from a unified view of the imaging data standards: DICOM, NIfTI, ANALYZE, AFNI, MINC, etc. The existence of a single package for handling imaging data formats would facilitate interoperability between the ever increasing number of \proglang{R} packages devoted to medical image analysis. We do not assume that the data structures in \pkg{oro.nifti} are best-suited for this purpose and we welcome an open discussion around how best to provide this standardization to the end user.
\section*{Acknowledgments}
The authors would like to thank the National Biomedical Imaging Archive (NBIA), the National Cancer Institute (NCI), the National Institute of Health (NIH) and all institutions that have contributed medical imaging data to the public domain. The authors would also like to thank K. Tabelow for providing functionality around the AFNI data format. VS is supported by the German Research Council (DFG SCHM 2747/1-1).
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.