\documentclass[article,nojss,shortnames]{jss}
%\VignetteIndexEntry{Using oro.nifti}
%\VignettePackage{NIfTI+ANALYZE+AFNI}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath,rotating}
\usepackage{thumbpdf}

%% almost as usual
\author{Brandon Whitcher\\Pfizer Worldwide R{\&}D \And 
        Volker J. Schmid\\Ludwig-Maximilians Universit\"at M\"unchen \AND
        Andrew Thornton\\Cardiff University}
\title{Working with the {NIfTI} Data Standard in \proglang{R}}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Brandon Whitcher, Volker J. Schmid, Andrew Thornton} %% comma-separated
\Plaintitle{Working with the {NIfTI} Data Standard in R} %% without formatting
\Shorttitle{{NIfTI} and R} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{

The package \pkg{oro.nifti} facilitates the interaction with and
manipulation of medical imaging data that conform to the ANALYZE,
NIfTI and AFNI formats.  The \proglang{S}4 class framework is used to
develop basic ANALYZE and NIfTI classes, where NIfTI extensions may be
used to extend the fixed-byte NIfTI header.  One example of this, that
has been implemented, is an \proglang{XML}-based ``audit trail''
tracking the history of operations applied to a data set.  The
conversion from DICOM to ANALYZE/NIfTI is straightforward using the
capabilities of \pkg{oro.dicom}.  The \proglang{S}4 classes have been
developed to provide a user-friendly interface to the ANALYZE/NIfTI
data formats; allowing easy data input, data output, image processing
and visualization.

}
\Keywords{export, imaging, import, medical, visualization}
\Plainkeywords{export, imaging, import, medical, visualization} %% without formatting
%% at least one keyword must be supplied

%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{13}
%% \Issue{9}
%% \Month{September}
%% \Year{2004}
%% \Submitdate{2004-09-29}
%% \Acceptdate{2004-09-29}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Brandon Whitcher\\
  Pfizer Worldwide Research \& Development\\
  610 Main Street\\
  Cambridge, MA 02139, United States\\
  E-mail: \email{bwhitcher@gmail.com}\\
  URL: \url{http://www.imperial.ac.uk/people/b.whitcher}, \url{http://rigorousanalytics.blogspot.com}\\
  
  Volker J. Schmid\\
  Bioimaging group\\
  Department of Statistics\\
  Ludwig-Maximilians-Universit\"at M\"unchen\\
  80539 M\"unchen, Germany\\
  E-mail: \email{volker.schmid@lmu.de}\\
  URL: \url{http://volkerschmid.de}\\
  
  Andrew Thornton\\
  Cardiff University School of Medicine\\
  Heath Park\\
  Cardiff CF14 4XN, United Kingdom\\
  E-mail: \email{art27@cantab.net}\\
  %% URL: \url{http://}\\
}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/1/31336-5053
%% Fax: +43/1/31336-734

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<preliminaries,echo=FALSE,results=hide>>=
library("bitops")
library("XML")
library("splines")
library("oro.nifti")
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
options(niftiAuditTrail = TRUE)
@

\begin{document}

%% include your article here, just as usual
%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.

\section{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(<obj>, "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 otuput.\\
      \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

<<mniLR_nifti>>=
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

<<mniLR-png,echo=FALSE,results=hide>>=
jpeg(filename="mniLR.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<mniLR-image>>=
image(mniLR)
@ 
<<mniLR-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

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]{mniLR.jpeg} &
      \includegraphics*[width=0.45\textwidth]{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

<<mniRL-read>>=
fname <- system.file(file.path("nifti", "mniRL.nii.gz"), package="oro.nifti")
(mniRL <- readNIfTI(fname))
@ 
<<mniRL-png,echo=FALSE,results=hide>>=
jpeg(filename="mniRL.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<mniRL-image>>=
image(mniRL)
<<mniRL-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

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]{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.

<<mniRL-ortho-png,echo=FALSE,results=hide>>=
jpeg(filename="mniRL_orthographic.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<mniRL-orthographic>>=
orthographic(mniRL)
@ 
<<mniRL-ortho-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

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

<<NIfTI-slots>>=
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.
<<NIfTI-constructor>>=
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).

<<NIfTI-write>>=
writeNIfTI(random.image, "random")
list.files(pattern="random")
@ 

\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

<<niftiAuditTrail,eval=FALSE>>=
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}}
<<NIfTI audit.trail 01>>=
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

<<EBImage01,eval=FALSE>>=
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]{ffd.jpeg}\\
      \textbf{(a)}\\
      \includegraphics*[width=0.6\textwidth]{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

<<ffd>>=
filtered.func.data <- 
  system.file(file.path("nifti", "filtered_func_data.nii.gz"), 
              package="oro.nifti")
(ffd <- readNIfTI(filtered.func.data))
@ 
<<ffd-png,echo=FALSE,results=hide>>=
jpeg(filename="ffd.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<ffd-image>>=
image(ffd, zlim=range(ffd)*0.95)
@ 
<<ffd-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

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.

<<ffd-ortho-png,echo=FALSE,results=hide>>=
jpeg(filename="ffd_orthographic.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<ffd-orthographic>>=
orthographic(ffd, xyz=c(34,29,10), zlim=range(ffd)*0.9)
@ 
<<ffd-ortho-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

\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]{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.

<<ffd-glm-design,echo=TRUE,results=hide>>=
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]
@ 
<<ffd-design.png,echo=FALSE,results=hide>>=
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.

<<ffd-glm,echo=TRUE,results=hide>>=
##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]{ffd_zstat1.jpeg}\\
      \textbf{(a)}\\
      \includegraphics*[width=0.6\textwidth]{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.

<<zstat1>>=
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)
@ 
<<zstat1-png,echo=FALSE,results=hide>>=
jpeg("ffd_zstat1.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<zstat1-overlay>>=
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)
<<zstat1-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 
<<zstat2-png,echo=FALSE,results=hide>>=
jpeg("ffd_zstat2.jpeg", width=480, height=480, quality=95, bg="black")
@ 
<<zstat2-overlay>>=
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)
<<zstat2-dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

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).

\bibliography{nifti}

\end{document}

