* Supposing another matrix of direction cosines is added to the itk::Image to represent the orientation of the measurement frmae, one still has to chose which frame is being used as the basis for the cosines- it is image frame or the world frame? Or should this be left as an application-specific choice?

* Supposing another matrix of direction cosines is added to the itk::Image to represent the orientation of the measurement frmae, one still has to chose which frame is being used as the basis for the cosines- it is image frame or the world frame? Or should this be left as an application-specific choice?

* When an image is rotated, who is responsible for updating vector- or tensor-valued pixels? Should the pixel values actually be changed, or is it enough to update a matrix representing the measurement frame?

* When an image is rotated, who is responsible for updating vector- or tensor-valued pixels? Should the pixel values actually be changed, or is it enough to update a matrix representing the measurement frame?

+

+

= Recursive Templates for fast computation of Index-Point transform =

+

+

Conversions betweeen index and physical point coordinates are critical for the performance of many ITK algorithms. It is desirable to speed up those operations even at the cost of some code complexity. The following is a proposal for a helper class that will compute the transformations between indices and physical points by taking advantage of the capabilities of C++ templates for performing recursion.

Proposal

itkImage should contain information regarding the orienttaion of the image volume.

Coordinate Systems

Coordinate systems are an important part of any medical application. Medical scanners create regular, "rectangular" arrays of points and cells. The topology of the arrays is implicit in the representation. The geometric location of each point is also implicit.

vtk and itk Image

vtk and itk store meta-information that can be used to convert between image coordinates (i-j-k) and world coordinates (x-y-z). Each image has an origin that is a 3-D (n-D for itk) vector. The origin (x-y-z) specifies the world coordinate of the first point in memory. The spacing specifies the distance between points along each axis. Using the spacing and origin, the transformation between i-j-k and x-y-z is a fast, simple computation.
x(yz) = origin + i(jk) * spacing;
i(jk) = (x(yz) - origin) / spacing;

Limitations of the current Image

Many image processing and segmentation algrothms do not need additional spatial information. But, registration and modeling techniques need to respect the orientation of the image arrays.

Proposed Extension to itk::Image

By adding an itk::Matrix containing direction cosines, the i-j-k -> x-y-z transformation could include the orientation in the computation. Adding the matrix will not change the existing API. All index to point calculations are confined to itk::Image. The transformations in matrix form are:
XYZ = To * Rc * Ss * IJK
where To is a translate to the origin, Rc is the matrix of direction cosines and Ss is a scale matrix of spacings.
There are performance considerations, so the implementation may cache some internal matrices and state.

Questions Raised by Proposal

The coordinate frame in which the direction cosines are measured needs to be either:

identified explicitly in the Image, so that the image itself knows if it is, for example, RAS (from NIFTY-1 data) versus LPS (from DICOM data)

not explicitly identified in the Image, but convertible from whatever is native to the input file, to whatever the application prefers to use.

Consensus on the 4 March 2005 tcon seemed to favor the second choice, using an (optional) orientation filter to impose one chosen orientation. This works most cleanly for axis aligned orientations, this is most commonly done for 3-D (4-D is also quite likely in fMRI, cardiac and other time based imaging methods).

http://wideman-one.com/gw/brain/analyze/formatdoc.htm
The ITK Analyze filter places itk::SpatialOrientation tags in the MetaDataDictionary. A utility class itkOrientImageFilter takes inputs of one 3D image, it's orientation, and the desired orientation and internally uses the flip and permute filters to acheive the desired behavior.

Some notes on the DICOM convention and current ITK usage

Image Orientation (Patient) (0020,0037)
The direction cosines of the first row and the first column with respect to
the patient.

Image Position (Patient) (0020,0032)
The x, y, and z coordinates of the upper left hand corner (center of the first
voxel transmitted) of the image, in mm.

DICOM Plane Attribute Descriptions C.7.6.2.1.1 Image
Position And Image Orientation
The Image Position (0020,0032) specifies the x, y, and z coordinates of the
upper left hand corner of the image; it is the center of the first voxel
transmitted.

Image Orientation (0020,0037) specifies the direction cosines of the
first row and the first column with respect to the patient. These Attributes
shall be provide as a pair. Row value for the x, y, and z axes respectively
followed by the Column value for the x, y, and z axes respectively.

The direction of the axes is defined fully by the patient's orientation.
The x-axis is increasing to the left hand side of the patient.
The y-axis is increasing to the posterior side of the patient.
The z-axis is increasing toward the head of the patient.
The patient based coordinate system is a right handed system, i.e. the vector
cross product of a unit vector along the positive x-axis and a unit vector
along the positive y-axis is equal to a unit vector along the positive z-axis.

Note: If a patient lies parallel to the ground, face-up on the table,
with his feet-to-head direction same as the front-to-back direction of the
imaging equipment, the direction of the axes of this patient based coordinate
system and the equipment based coordinate system in previous versions of this
Standard will coincide. The Image Plane Attributes, in conjunction with the
Pixel Spacing Attribute, describe the position and orientation of the image
slices relative to the patient-based coordinate system.
In each image frame the Image Position (Patient) (0020,0032) specifies the
origin of the image with respect to the patient-based coordinate system.

Reference Coordinate System (RCS) and the Image Orientation (Patient)
(0020,0037) attribute values specify the orientation of the image frame rows
and columns. The mapping of pixel location (i,j) to the RCS is calculated as
follows:

2) The row and column direction cosine vectors shall be normal, i.e. the dot product of each direction cosine vector with itself shall be unity.

The direction cosine for slices (k) is determined by the cross product of the
two provided direction cosines. The coordinate system is right handed and the
direction cosine must be unit magnitude.

Determination of whether consecutive slices were acquired in a +Z or a -Z
direction requires comparison of the Image Position (Patient) tag values
of at least two slices.

Some other useful DICOM tags include:

Pixel Spacing (0028,0030) Physical distance in the patient between the
center of each pixel, specified by a numeric pair - adjacent row spacing
(delimiter) adjacent column spacing in mm.

Slice Thickness (0018,0050) Nominal slice thickness, in mm.
Note that the value of Slice Thickness can be equal to, larger than, or
smaller than, the distance between adjacent slices.

DICOM on Diffusion MRI

Diffusion Gradient Orientation (0018, 9089) - The direction cosines of the
diffusion gradient vector with respect to the patient. Required if Frame Type (0008,9007) Value 1 of this frame is ORIGINAL. May be present otherwise.

Current ITK Usage and sources of confusion

It seems unfortunate to have a letter code selection for axes
that is opposite to that used by DICOM and hence the scanner manufacturers.

itkSpatialOrientation.h -

// Coordinate orientation codes have a place-value organization such that
// an ImageDimension-al sequence of subcodes says both which varies
// fastest
// through which varies slowest, but also which end of the frame of
// reference
// is considered zero for each of the coordinates. For example, 'RIP'
// means
// Right to Left varies fastest, then Inferior to Superior, and Posterior
// to
// Anterior varies the slowest.

I take this to mean R implies the fastest varying axis runs from Right to
Left. DICOM would give this the code L, as described here:

C.7.6.1.1.1 Patient Orientation The Patient Orientation (0020,0020) relative to the image plane shall be specified by two values that designate the anatomical direction of the positive row axis (left to right) and the positive column axis (top to bottom). The first entry is the direction of the rows, given by the direction of the last pixel in the first row from the first pixel in that row. The second entry is the direction of the columns, given by the direction of the last pixel in the first column from the first pixel in that column. Anatomical direction shall be designated by the capital letters: A (anterior), P (posterior), R (right), L (left), H (head), F (foot). Each value of the orientation attribute shall contain at least one of these characters. If refinements in the orientation descriptions are to be specified, then they shall be designated by one or two additional letters in each value. Within each value, the letters shall be ordered with the principal orientation designated in the first character.

Interpretation in DICOM terms:
orientation is CORONAL, it turns out slices running from A to P.
unit vector along row is L i.e. a unit vector along the row points from patient Right to Patient Left.
unit vector along column is I (i.e. columns from the Superior to the Inferior)

row : L
columns : I
slices : P

Three letter code for (column index, row index, slice index is) ILP.
This indicates the origin is at the superior, right, anterior corner of the
volume, and that the axes run from superior to inferior, from right to left,
from anterior to posterior.

This is the OPPOSITE CONVENTION TO THAT DESCRIBED BY METAIO
AnatomicalOrientation tag. [This is being fixed in the documentation.
The implementation follows the DICOM convention - see itkMetaImageIO]

Note that data described as RAS using this convention is not the same as data
described as RAS using the DICOM convention.

Similarly, the InsightSNAP tool in InsightApplications (which has a great IO
loader by the way), uses the following:
When it says LPS it means the x axis is RUNNING FROM left to right,
the y axis is running from Posterior to Anterior, the z axis is running from
Superior to Inferior.

This is exactly the same as what DICOM means when it says RAI.

The third coordinate system: measurement or laboratory frame

Whenever non-scalar pixel values involve vectors (or tensors) with coefficients that are expressed with respect to some coordinate frame, the orientation of that coordinate frame should also be known. This frame, which can be called the measurement frame, is conceptually distinct from both the "world" space around the array (such as right-anterior-superior), and the image-aligned basis determined by the direction cosines.

With medical data, the measurement frame becomes an issue for expressing the diffusion-sensitizing gradient directions used in diffusion-weighted (DW) MRI: the diffusion tensors estimated from DW-MRI will likely be expressed in whatever frame the gradients were expressed in. This frame could be, for example, the same as the one used for world space, or the image space frame. The DICOM spec includes a field (see notes above) for identifying the space of the gradients, but apparently scanners are not consistently responsible in their setting of this field, so it may have to be manually recorded outside the DICOM file.

But outside the specific context of DW-MRI and DT-MRI, the measurement frame is an issues for any kind of vector and tensor data, so perhaps it should have representation/handling in ITK.

Are there already established conventions (from ITK users/code) for the frame in which ITK vectors are expressed?

Supposing another matrix of direction cosines is added to the itk::Image to represent the orientation of the measurement frmae, one still has to chose which frame is being used as the basis for the cosines- it is image frame or the world frame? Or should this be left as an application-specific choice?

When an image is rotated, who is responsible for updating vector- or tensor-valued pixels? Should the pixel values actually be changed, or is it enough to update a matrix representing the measurement frame?

Recursive Templates for fast computation of Index-Point transform

Conversions betweeen index and physical point coordinates are critical for the performance of many ITK algorithms. It is desirable to speed up those operations even at the cost of some code complexity. The following is a proposal for a helper class that will compute the transformations between indices and physical points by taking advantage of the capabilities of C++ templates for performing recursion.