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 gg/mri_g++-* library and its header is <gg/mri/image.H>. A number of 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 1 in the fourth dimension, while a DTI will have resolution 6 in the fourth dimension.
• i.n: a raw 2-D image slice, stored in fixed-point format as little-endian uint16. 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 printf("%03d",n); 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:
  • i.001 stores the slice for Z=0 and the first tensor component
  • i.002 stores the slice for Z=0 and the second tensor component
  • ...
  • i.006 stores the slice for Z=0 and the sixth tensor component
  • i.007 stores the slice for Z=1 and the first tensor component
  • ...
  • i.012 stores the slice for Z=1 and the sixth tensor component
  • ..
  • i.n stores the slice for $Z=\lfloor (n-1)/6\rfloor$ and the $((n-1)mod6)+1$-th tensor component
• vsize: a text file indicating the voxel size in millimeters
• parameters: a text file indicating the fixed-point parameters (scaleIntensity and offsetIntensity) and imaging parameters such as TR, TE, and data-to-lab-frame transformation

Memory Format

The uint16 data value at index (i,j,k,d) is stored at linear index i + j*xres + k*xres*yres + d*xres*yres*zres in the data array. The uint16 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 offset and offset + scale * 65535). 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:

• image_value = offset + scale * data_value
• data_value = (image_value - offset) / scale

The scale and offset values may be read with getScaleOffset() and set with setScaleOffset().

The functions for data access specify coordinates as (x,y,z,d), 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
  • void set(int x, int y, int z, int d, double val) --- set the data value at index (x,y,z,d) to val. Even though val is a double, 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
  • double lookup(int x, int y, int z, int d, double bg, int *inside) --- read the data value at index (x,y,z,d), returning bg if the indices are out of bounds, and setting inside to 0 in that case.
  • double lookup_s(int x, int y, int z, int d, double bg, int *inside) --- read the image value at index (x,y,z,d), returning bg if the indices are out of bounds, and setting inside to 0 in that case. The "_s" at the end of the function name stands for "scaled".
  • double lookup_nocheck(int x, int y, int z, int d) --- read the data value at index (x,y,z,d), with no bounds checking. You'll get memory errors (or random junk) if you read out of bounds.
  • double lookup_nocheck(int x, int y, int z, int d) --- read the image value at index (x,y,z,d), with no bounds checking. You'll get memory errors (or random junk) if you read out of bounds.
• Interpolation
  • double sample(interp, double x, double y, double z, int d, double back=0, double *inside=0, double *g=0) --- smoothly sample (that is, interpolate) the data value at fractional-valued index (x,y,z,d) (d may not be fractional-valued; no interpolation occurs along the fourth dimension), assuming that out-of-bounds voxels have value back. Non-null values of g are not supported, so just stick with the default (I'm not even really sure what that argument is supposed to do). The interp argument is the interpolation kernel keyword, the options for which are:
    • LINEAR --- linear interpolation among the eight corners of the containing cube
    • CUBIC --- cubic spline interpolation
    • PFIT --- quadratic interpolation
    • CONST_INTERP --- nearest-neighbor
    • CONST_NOCHECK --- nearest-neighbor, with no bounds-checking
  • double sample(interp, double x, double y, double z, int d, double back=0, double *inside=0, double *g=0) --- 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 (i,j,k) to spatial coordinates (x,y,z) 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: x = i * vsize[0] + space_shift_x; y = j * vsize[1] + space_shift_y; z = k * vsize[3] + space_shift_z; (That is, data indices (i,j,k) map directly to spatial coordinates (x,y,z), 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:
  • +x = Patient's Left (x increases from patient's right to patient's left)
  • +y = Posterior (y increases from anterior (nose) to posterior (back of head))
  • +z = Superior (z increases from inferior (base of head) to superior (top of head))

Technically, the setXform() and getXform() 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.