View source for MRIimage

MRIimage is our in-house format for storing MRI data and is used as the underlying file format for most of what we do in the [[Diffusion MRI]] project.  It is manipulated by the <code>gg/mri_g++-*</code> library and its header is <code>&lt;gg/mri/image.H&gt;</code>.  A number of [[diffusion MRI Matlab scripts|Matlab scripts]] have also been written to interact with this file format.

== File Format ==

In MRIimage, an image volume is represented by a directory containing a number of files:
* '''resolution''': a text file indicating the resolution of the image in voxels.  Images are 4-D; a single DWI will have resolution <tt>1</tt> in the fourth dimension, while a DTI will have resolution <tt>6</tt> in the fourth dimension.
* '''i.''n''''': a raw 2-D image slice, stored in fixed-point format as little-endian <tt>uint16</tt>.  ''n'' corresponds to the X-Y slice number, which increases sequentially as the index in the fourth dimension increases, and is formatted as a string in the filename like <code>printf("%03d",n)</code>; hence numbers <1000 are zero-padded to three digits.  Since the fourth dimension is of unit size in a DWI, ''n'' increases sequentially as the index in Z increases.  For a DTI, however:
** <tt>i.001</tt> stores the slice for Z=0 and the first tensor component
** <tt>i.002</tt> stores the slice for Z=0 and the second tensor component
** ...
** <tt>i.006</tt> stores the slice for Z=0 and the sixth tensor component
** <tt>i.007</tt> stores the slice for Z=1 and the first tensor component
** ...
** <tt>i.012</tt> stores the slice for Z=1 and the sixth tensor component
** ..
** <tt>i.''n''</tt> stores the slice for <math>Z = \lfloor (n-1)/6 \rfloor</math> and the <math>((n-1) \mod 6) + 1</math>-th tensor component
* '''vsize''': a text file indicating the voxel size in millimeters
* '''parameters''': a text file indicating the fixed-point parameters (<tt>scaleIntensity</tt> and <tt>offsetIntensity</tt>) and imaging parameters such as TR, TE, and data-to-lab-frame transformation

== Memory Format ==
The <tt>uint16</tt> data value at index <tt>(i,j,k,d)</tt> is stored at linear index <code>i + j*xres + k*xres*yres + d*xres*yres*zres</code> in the data array.  The <tt>uint16</tt> value stored at this location (which is an integer between 0 and 65535) represents an actual image value in a fixed-point format (which may be a fractional-valued scalar between <tt>offset</tt> and <tt>offset + scale * 65535</tt>).  In the discussion below, values of the former type are called "data values", while the latter type are called "image values".  The conversions between data values and image values are as follows:
* <code>image_value = offset + scale * data_value</code>
* <code>data_value = (image_value - offset) / scale</code>
The scale and offset values may be read with <code>getScaleOffset()</code> and set with <code>setScaleOffset()</code>.

The functions for data access specify coordinates as <tt>(x,y,z,d)</tt>, but these values actually refer to voxel indices, not spatial coordinates, which are described in the next section.

The data access functions in the MRIimage library are as follows:
* '''Writing'''
** <code>void set(int x, int y, int z, int d, double val)</code> --- set the ''data'' value at index <tt>(x,y,z,d)</tt> to <tt>val</tt>.  Even though <tt>val</tt> is a <tt>double</tt>, it is still a ''data'' value; it will be rounded to the nearest integer and truncated to [0, 65535] before writing to the data array.
* '''Reading'''
** <code>double lookup(int x, int y, int z, int d, double bg, int *inside)</code> --- read the ''data'' value at index <tt>(x,y,z,d)</tt>, returning <tt>bg</tt> if the indices are out of bounds, and setting <tt>inside</tt> to <tt>0</tt> in that case.
** <code>double lookup_s(int x, int y, int z, int d, double bg, int *inside)</code> --- read the ''image'' value at index <tt>(x,y,z,d)</tt>, returning <tt>bg</tt> if the indices are out of bounds, and setting <tt>inside</tt> to <tt>0</tt> in that case.  The "<tt>_s</tt>" at the end of the function name stands for "scaled".
** <code>double lookup_nocheck(int x, int y, int z, int d)</code> --- read the ''data'' value at index <tt>(x,y,z,d)</tt>, with no bounds checking.  You'll get memory errors (or random junk) if you read out of bounds.
** <code>double lookup_nocheck(int x, int y, int z, int d)</code> --- read the ''image'' value at index <tt>(x,y,z,d)</tt>, with no bounds checking.  You'll get memory errors (or random junk) if you read out of bounds.
* '''Interpolation'''
** <code>double sample(interp, double x, double y, double z, int d, double back=0, double *inside=0, double *g=0)</code> --- smoothly sample (that is, interpolate) the ''data'' value at fractional-valued index <tt>(x,y,z,d)</tt> (<tt>d</tt> may not be fractional-valued; no interpolation occurs along the fourth dimension), assuming that out-of-bounds voxels have value <tt>back</tt>.  Non-null values of <tt>g</tt> are not supported, so just stick with the default (I'm not even really sure what that argument is supposed to do).  The <tt>interp</tt> argument is the interpolation kernel keyword, the options for which are:
*** <tt>LINEAR</tt> --- linear interpolation among the eight corners of the containing cube
*** <tt>CUBIC</tt> --- cubic spline interpolation
*** <tt>PFIT</tt> --- quadratic interpolation
*** <tt>CONST_INTERP</tt> --- nearest-neighbor
*** <tt>CONST_NOCHECK</tt> --- nearest-neighbor, with no bounds-checking
** <code>double sample(interp, double x, double y, double z, int d, double back=0, double *inside=0, double *g=0)</code> --- same as above, but return the ''image'' value rather than the ''data'' value.

== Coordinate System ==
For any image format that represents real-world measurements, the mapping from data array indices <tt>(i,j,k)</tt> to spatial coordinates <tt>(x,y,z)</tt> must be defined.  The axes of this spatial coordinate system must also be defined.  More mature formats like [[NIfTI]] are explicit about both this mapping and the axes, but MRIimage works by convention.  The convention for the MRIimage library (that is, the representation of the image in memory) is this:

* '''Mapping''': <code>x = i * vsize[0] + space_shift_x;</code> <code>y = j * vsize[1] + space_shift_y;</code> <code>z = k * vsize[3] + space_shift_z;</code> (That is, data indices <tt>(i,j,k)</tt> map directly to spatial coordinates <tt>(x,y,z)</tt>, with voxel sizes and the space shift factored in.)
* '''Axes''': The axes of the spatial coordinate system are a compromise between conventions for clinical imaging and computer graphics.  The axes are defined as:
** <tt>+x</tt> = Patient's Left (<tt>x</tt> increases from patient's right to patient's left)
** <tt>+y</tt> = Posterior (<tt>y</tt> increases from anterior (nose) to posterior (back of head))
** <tt>+z</tt> = Superior (<tt>z</tt> increases from inferior (base of head) to superior (top of head))

Technically, the <tt>setXform()</tt> and <tt>getXform()</tt> functions may be used to change the mapping from data indices to spatial coordinates from the default vsize- and space_shift-based mapping, but none of the code in the [[diffusion processing pipeline]] uses this feature, and it may not be fully implemented in all programs.

[[Category:Diffusion MRI]][[Category:File Formats]]