Working with the DICOM and NIfTI Data Standards in R

Two packages, oro.dicom and oro.nifti, are provided for the interaction with and manipulation of medical imaging data that conform to the DICOM standard or ANALYZE/NIfTI formats. DICOM data, from a single file or directory tree, may be uploaded into R using basic data structures: a data frame for the header information and a matrix for the image data. A list structure is used to organize multiple DICOM files. The S4 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 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 both packages. The S4 classes have been developed to provide a userfriendly interface to the ANALYZE/NIfTI data formats; allowing easy data input, data output, image processing and visualization.


Introduction
Medical imaging is well established in both the clinical and research areas with numerous equipment manufacturers supplying a wide variety of modalities. The DICOM (Digital Imaging and Communications in Medicine; http://medical.nema.org) standard was developed from earlier standards and released in 1993. It is the data format for clinical imaging equipment and a variety of other devices whose complete specification is beyond the scope of this paper. All major manufacturers of medical imaging equipment (e.g., GE, Siemens, Philips) have so-called DICOM conformance statements that explicitly state how their hardware implements DICOM. The DICOM standard provides interoperability across hardware, but was not designed to facilitate efficient data manipulation and image processing. Hence, additional data formats have been developed over the years to accommodate data analysis and image processing.
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.
The material presented here provides users with a method of interacting with DICOM, ANA-LYZE and NIfTI files in R (R Development Core Team 2010). Real-world data sets, that are publicly available, are used to illustrate the basic functionality of the two packages: oro.dicom and oro.nifti. It should be noted that both packages focus on functions for data input/output and visualization. S4 classes nifti and anlz are provided for further statistical analysis in 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 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 nifti and anlz objects. Additionally, the oro.nifti package allows one to track every operation on a nifti object in an XML-based audit trail.
The oro.dicom and oro.nifti packages should appeal not only to R package developers, but also to scientists and researchers who want to interrogate medical imaging data using the statistical capabilities of R without writing and validating their own basic data input/output functionality. Table 1 lists the key functions for both packages and groups them according to common functionality. An example of using statistical methodology in R for the analysis of functional magnetic resonance imaging (fMRI) data is given in section 3.7. Packages already available on CRAN that utilize oro.dicom and oro.nifti include cudaBayesreg (Ferreira da Silva 2010a), dcemriS4 (Whitcher and Schmid 2011), dpmixsim (Ferreira da Silva 2010b), and RNiftyReg (Clayden 2011).
2. oro.dicom: DICOM data input/output in R The DICOM "standard" for data acquired using a clinical imaging device is very broad and complex. Roughly speaking each DICOM-compliant file is a collection of fields organized into two two-byte sequences (group,element) Data element with explicit VR other than as shown above: Data element with implicit VR: Figure 1: Byte ordering for a single (group,element) tag in the DICOM standard. Explicit VRs store the VR as text characters in two bytes. More information is provided in Section 7, Part 3.5-2009 of the DICOM standard (http://medical.nema.org).
forthcoming in the file. There is no fixed number of bytes for a DICOM header. The final (group,element) tag should be the "pixel data" tag (7FE0,0010), such that all subsequent information is related to the image(s).
All attributes in the DICOM standard require different data types for correct representation. These are known as value representations (VRs) in DICOM, which may be encoded explicitly or implicitly. There are 27 explicit VRs defined in the DICOM standard. Detailed explanations of these data types are provided in the Section 6.2 (part 5) of the DICOM standard (http://medical.nema.org). Internal functions have been written to manipulate each of the value representations and are beyond the scope of this article. The functions str2date and str2time are useful for converting from the DICOM Datetime and Time value representations to R date and time objects, respectively.

The DICOM header
Accessing the information stored in a single DICOM file is provided using the dicomInfo function. The basic structure of a DICOM file is summarized in Figure 1, for both explicit and implicit value representations. The first two bytes represent the group tag and the second two bytes represent the element tag, regardless of the type of VR. The third set of two bytes contains the characters of the VR on which a decision about being implicit or explicit is made. Explicit VRs of type (OB, OF, OW, SQ, UT, UN) skip bytes six and seven (counting from zero), convert the next four bytes into an integer length and read length number of objects from the DICOM file. All other explicit VRs follow a slightly different path where bytes six and seven (counting from zero) provide an integer length and all remaining bytes are read in as the value. If the character string in bytes four and five do not correspond to a known VR (Figure 1), then the (group,element) tag is declared to be implicit, the length is taken from bytes four through seven and all remaining bytes contribute to the value.
The basic structure of the resulting object is a list with two elements: the DICOM header (hdr) and the DICOM image (img). The header information is organized in a data frame with six columns and an unknown number of rows depending on the input parameters.
R> fname <-system.file(file.path("dcm", "Abdo.dcm"), package = "oro.dicom") R> abdo <-dicomInfo(fname) R> names(abdo) [1] "hdr" "img" The first five columns are taken directly from the DICOM header information (group, element, code, length and value) or inferred from that information (name). Note, the (group,element) values are stored as character strings even though they are hexadecimal numbers. All aspects of the data frame may be interrogated in R in order to extract relevant information from the DICOM header; e.g., "BitsAllocated" as above. The sequence column is used to keep track of tags that are embedded in a fixed-length SequenceItems tag or between a SequenceItem-SequenceDelimitationItem pair. When multiple DICOM files are located in a single directory, or spread across multiple directories, one may use the function dicomSeparate (applied here to the directory hk-40).

R> unlist(lapply(hk40, length))
hdr img 40 40 The object associated with dicomSeparate is now a nested set of lists, where the hdr element is a list of data frames and the img element is a list of matrices. These two lists are associated in a pairwise sense; i.e., hdr[ [1]] is the header information for the image img[ [1]]. Default parameters recursive = TRUE and pixelData = TRUE (which is actually an input parameter for dicomInfo) allow the user to search down all possible sub-directories and upload the image in addition to the header information, respectively. Also, by default all files are treated as DICOM files unless the exclude parameter is set to the unwanted file extension; e.g., exclude = "xml".
The list of DICOM header information across multiple files may be converted to a single data frame using dicomTable, and written to disc for further analysis; e.g., using write.csv.
R> unique(extractHeader(hk40$hdr, "SliceThickness")) The tag SliceLocation is extracted from the DICOM header information (at the first element in the list) and processed using the diff function, and should agree with the SliceThickness tag. Single DICOM fields may also be extracted from the list of DICOM header information that contain attributes that are crucial for further image processing; e.g., extracting relevant MR sequences or acquisition timings.

The DICOM image
Most DICOM files involve a single slice from an acquisition -the image. A notable exception is the Siemens MOSAIC format (addressed in Section 2.2.1). The oro.dicom package assumes the image is stored as a flat file of two-byte integers without compression. A variety of additional image formats are possible within the DICOM standard; e.g., RGB-colorized, JPEG, JPEG Lossless, JPEG 2000 and run-length encoding (RLE). None of these formats are currently available in oro.dicom. Going back to the Abdo.dcm example, the image is accessed via R> image(t(abdo$img), col = grey(0:64/64), axes = FALSE, xlab = "", + ylab = "") where the transpose operation is necessary for proper visualization of the image. Figure 2 displays a coronal slice through the abdomen from an MRI acquisition. All information from the original data acquisition should accompany the image through the DICOM header, and this information is utilized as much as possible by oro.dicom to simplify the manipulation of DICOM data. As previously shown, this information is easily available to the user by matching DICOM header fields with valid strings. Note, the function extractHeader assumes the output should be coerced via as.numeric but this may be disabled setting the input parameter numeric=FALSE.
R> extractHeader(abdo$hdr, "Manufacturer", numeric = FALSE) R> extractHeader(abdo$hdr, "RepetitionTime") R> extractHeader(abdo$hdr, "EchoTime") [1] 100 The basic DICOM file structure does not encourage the analysis of multi-dimensional imaging data (e.g., 3D or 4D) commonly acquired on clinical scanners. Hence, the oro.dicom package has been developed to access DICOM files and facilitate their conversion to the NIfTI or ANALYZE formats in R. The conversion process requires the oro.nifti package and will be outlined in Section 4.

Siemens MOSAIC format
Siemens multi-slice EPI (echo planar imaging) data may be collected as a "mosaic" image; i.e., all slices acquired in a single TR (repetition time) frame of a dynamic run are stored in a single DICOM file. The images are stored in an M ×N array of images. The function create3D will try to guess the number of images embedded within the single DICOM file using the AcquisitionMatrix field. If this doesn't work, one may enter the (M, N ) doublet explicitly.
R> fname <-system.file(file.path("dcm", "MR-sonata-3D-as-Tile.dcm"),   Figure 3b is displayed after re-organizing the original DICOM file into a threedimensional array (it was also converted to the NIfTI format for ease of visualization using the overloaded image function in oro.nifti).
3. oro.nifti: NIfTI-1 data input/output in 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 oro.nifti (accessed via system.file("doc/ANALYZE75.pdf", package="oro.dicom")) since it does not appear to be available from 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; http://nifti.nimh.nih.gov). The NIfTI-1 data format is almost identical to the ANALYZE format, but offers a few improvements merging of the header and image information into one file (.nii) re-organization of the 348-byte fixed header into more relevant categories possibility of extending the header information.

The 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 (magic = "n+1\0") or two separate files (magic = "ni1\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 avg152T1_LR_nifti.nii.gz was obtained from the NIfTI website (http://nifti.nimh. nih.gov/nifti-1/). Successful execution of the commands R> fname <-system.file(file.path("nifti", "mniLR.nii.gz"), + package = "oro.nifti") R> (mniLR <-readNIfTI(fname)) R> descrip(mniLR) [1] "FSL3.2beta" produces an S4 "nifti" object (or "niftiAuditTrail" if the audit trail option is set). Two accessor functions are also provided: aux.file and 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).

The NIfTI image
Image information begins at the byte position determined by the voxoffset slot. In a single NIfTI file (magic = "n+1\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 voxoffset for single NIfTI files. The split NIfTI (magic = "ni1\0") and ANA-LYZE formats contain pairs of files, where the header and image information are separated, and do not have this problem. In this case voxoffset is set to 0.
The image function has been overloaded so that it behaves differently when dealing with medical image objects (nifti and anlz). The command

R> image(mniLR)
produces a three-dimensional array of the MNI brain, with the default NIfTI axes, and is displayed on a 10×10 grid of images ( Figure 4a). The image function for medical image S4 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 plot.type="single" in conjunction with the option z= to specify the slice. The second example of reading in and displaying medical imaging data in the NIfTI format avg152T1_RL_nifti.nii.gz was also obtained from the NIfTI website (http://nifti.nimh. nih.gov/nifti-1/). Successful execution of the commands R> fname <-system.file(file.path("nifti", "mniRL.nii.gz"), + package = "oro.nifti") R> (mniRL <-readNIfTI(fname))

NIfTI-1 format
Type : niftiAuditTrail Data Type : 2 (UINT8)  produces a three-dimensional array of the MNI brain that is displayed in a 10×10 grid of images ( Figure 4b). The two sets of data in Figure 4 are stored in two different orientations, commonly referred to as the radiological and 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.
An additional graphical display function has been added for nifti and anlz objects that allows a so-called orthographic visualization of the data.

R> orthographic(mniRL)
As seen in Figure 5 the mid-axial, mid-sagittal and mid-coronal planes are displayed by default. The slices used may be set using 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.

A note on axes and orientation
The NIfTI format contains an implicit generalized spatial transformation from the data coordinate system (i, j, k) into a real-space "right-handed" co-ordinate system. In this real-space system, the (x, y, z) axes are 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 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 oro.nifti will attempt to address more complicated spatial transformations and provide functionality to display the (x, y, z) axes on orthographic plots.
As introduced in Section 3.1 there are currently only two accessor functions to slots in the NIfTI header (aux.file and 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 nifti produces valid NIfTI objects, including a consistent header, from an arbitrary array.
R> n <-100 R> (random.image <-nifti(array(runif(n * n), c(n, n, 1)))) The function writeNIfTI outputs valid NIfTI class files, which can be opened in other medical imaging software. Files can either be stored as standard .nii files or compressed with gnuzip (default).

The audit trail
Following on from the S4 implementation of both the NIfTI and ANALYZE data formats, the ability to extend the NIfTI data format header is utilized in the oro.nifti package. First, extensions are properly handled when reading and writing NIfTI data. Second, users are allowed to add extensions to newly-created NIfTI objects by casting them as niftiExtension objects and adding niftiExtensionSection objects to the extensions slot. Third, by default all operations that are performed on a NIfTI object will generate what we call an audit trail that consists of an XML-based log (Temple Lang 2010). 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 R, user name, date, time, etc. When writing NIfTI-class objects to disk, the XML-based NIfTI extension is converted into plain text and saved using ecode=6 to denote plain ASCII text only. The user may control the tracking of data manipulation via the audit trail using a global option. For example, please use the command

R> options(niftiAuditTrail = FALSE)
to turn off the "audit trail" option in oro.nifti. Figure 6 displays output from the accessor function audit.trail(mniLR), the XML-based audit trail that is stored as a NIfTI header extension.

Interactive visualization
Basic visualization of nifti and anlz class images can be achieved with any visualization for arrays in R. For example, the EBImage package provides functions display and animate for visualization (Sklyar et al. 2010). Please note that functions in EBImage expect grey-scale values in the range [0, 1], hence the display of nifti data may be performed using R> mniLR.range <-range(mniLR) R> 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 R at this point in time. The mritc package provides basic interactive visualization of ANALYZE/NIfTI data using a Tcl/Tk interface (Feng and Tierney 2010).

An example using functional MRI data
This is an example of reading in, and displaying, a four-dimensional medical imaging data set in NIfTI format filtered_func_data obtained from the FSL evaluation and example data suite (http://www.fmrib.ox.ac.uk/fsl/fsl/feeds.html). Successful execution of the commands R> (ffd <-readNIfTI("filtered_func_data")) produces a four-dimensional (4D) array of imaging data that may be displayed in a 5×5 grid of images (Figure 7a). 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×4×6 mm) and reasonable temporal resolution (3 s). Figure 7b depicts the orthographic display of the filtered_func_data using the axial plane containing the left-and-right thalamus to approximately center the crosshair vertically.

Statistical analysis
The 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 (Friston et al. 1994(Friston et al. , 1995. The 4D volume of imaging data in filtered_func_data was acquired in an experiment with a repetition time TR = 3 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 µ = 6 and standard deviation σ = 3, is utilized here which is similar to the default values in FSL (Smith et al. 2004). 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. Figure 8 depicts the visual and auditory stimuli, convolved with the HRF, in the order of acquisition. The design matrix is than used in a voxel-wise GLM, where the lsfit function in 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 ls.print function. 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 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.
R> dof <-ntim(ffd) -1 R> Z.visual <-nifti(qnorm(pt(ffd.glm [1, , , ], dof, log.p = TRUE), + log.p = TRUE), datatype = 16) R> Z.auditory <-nifti(qnorm(pt(ffd.glm[2, , , ] 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 overlay command in 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 9. 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 9a) and the primary auditory cortex in the temporal lobe shows activation under the auditory stimulus (Figure 9b).

Converting DICOM to NIfTI
The oro.dicom and oro.nifti packages have been specifically designed to use as much information as possible from the metadata-rich DICOM format and apply that information in the construction of the NIfTI data volume. The function dicom2nifti converts a list of DICOM images into an nifti object, and likewise dicom2analyze converts such a list into an anlz object.
Historically, data conversion from DICOM to NIfTI (or ANALYZE) has been provided outside of R using one of several standalone software packages: Xmedcon (Nolf 2003),

An example using a single-series data set
Using the 40 images from the hk40 object (previously defined in Section 2.1) it is straightforward to perform DICOM-to-NIfTI conversion using only default settings and plot the results in either lightbox or orthographic displays.
The functions dicom2nifti and dicom2analyze will fail when the dimensions of the individual images in the DICOM list do not match. However, they do not check for different series numbers or patient IDs so caution should be exercised when scripting automated workflows for DICOM-to-NIfTI conversion. In cases where a DICOM file includes images from more than one series, the corresponding slices have to be chosen before conversion, using dicomTable, extractHeader, and matchHeader.

An example using a multiple-volume data set
The National Biomedical Imaging Archive (NBIA; http://cabig.nci.nih.gov/tools/NCIA) is a searchable, national repository integrating in vivo cancer images with clinical and genomic data. The NBIA provides the scientific community with public access to DICOM images, image markup, annotations, and rich metadata. The multiple MRI sequences processed here were downloaded from the "RIDER Neuro MRI" collection at http://wiki.nci.nih.gov/ display/CIP/RIDER. A small for loop has been written to operate on a subset of the DI-COM directory structure, where the SeriesInstanceUID DICOM header field is assumed to be 100% accurate in series differentiation.
There is always a balance between what information should be pre-specified versus what can easily be extracted from the DICOM headers or images.

Conclusion
Medical image analysis depends on the efficient manipulation and conversion of DICOM data. The oro.dicom and oro.nifti packages have 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 R will benefit from a unified view of the imaging data standards: DICOM, NIfTI and ANALYZE. The existence of a single package for handling imaging data formats would facilitate interoperability between the ever increasing number of R packages devoted to medical image analysis. We do not assume that the data structures in oro.dicom or 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.