Proposals:Orientation
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).
External References
- Current ITK SpatialOrientation header:
- How NRRD handles orientation information:
http://teem.sourceforge.net/nrrd/format.html#space
- Some info on Analyze format:
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.
/**
* \enum ValidAnalyzeOrientationFlags
* Valid Orientation values for objects
* - Key Description Origin dims[1] dims[2] dims[3]
* - =================================================================
* - 0 transverse-unflipped IRP R->L P->A I->S
* - 1 coronal-unflipped IRP R->L I->S P->A
* - 2 sagittal-unflipped IRP P->A I->S R->L
* - 3 transverse-flipped IRA R->L A->P I->S
* - 4 coronal-flipped SRP R->L S->I P->A
* - 5 sagittal-flipped ILP P->A I->S L->R
* - Where the Origin disignators are with respect to the patient
* - [(I)nferior|(S)uperior] [(L}eft|(R)ight] [(A)nterior|(P)osterior]
* \note Key's 0-5 correspond to the Analyze v7.5 orientations, and should not be changed.
*/
switch (temporient)
{
case itk::AnalyzeImageIO::ITK_ANALYZE_ORIENTATION_RPI_TRANSVERSE:
coord_orient = itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI;
}
itk::EncapsulateMetaData<itk::SpatialOrientation::ValidCoordinateOrientationFlags>(thisDic,ITK_CoordinateOrientation, coord_orient);
Code/Common/itkSpatialOrientation.h
Code/BasicFilters/itkOrientImageFilter.h
orienter->SetDesiredCoordinateOrientation( DesiredImageOrientation ); orienter->SetGivenCoordinateOrientation( fixedOrientation ); orienter->SetInput( inputFixedImage );
- DICOM can be deciphered in a way similar to the Analyze method as described by David Clunie in this mail message
- Info from Tosa Yasunari about how Freesurfer handles coordinates:
http://www.nmr.mgh.harvard.edu/~tosa/#coords
- A lot of good documentation is included in the nifti1.h header file:
http://nifti.nimh.nih.gov/pub/dist/src/nifti1.h
- The DICOM Standard in PDF:
http://medical.nema.org/dicom/2004.html Look in particular at http://medical.nema.org/dicom/2004/04_04PU.PDF section C.7.6.2.1.1, page 275 for a description of the coordinate transformation details.
Some notes on the DICOM convention and current ITK usage
- DICOM
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:
P_x = [ X_x \delta_i Y_x \delta_j 0 S_x ] [i ]
P_y [ X_y \delta_i Y_y \delta_j 0 S_y ] [j ]
P_z [ X_z \delta_i Y_z \delta_j 0 S_z ] [0 ]
1 [ 0 0 0 1 ] [1 ]
where :
Pxyz - The coordinates of the voxel (i,j) in the frame's image plane in units of mm.
Sxyz - The three values of the Image Position (Patient) (0020,0032) attributes. It is the location in mm from the origin of the RCS.
Xxyz - The values from the row (X) direction cosine of the Image Orientation (Patient) (0020,0037) attribute.
Yxyz - The values from the column (Y) direction cosine of the Image Orientation (Patient) (0020,0037) attribute.
i - Column index to the image plane. The first column is index zero.
\delta_i - Column pixel resolution of the Pixel Spacing (0028,0030) attribute in units of mm.
j - Row index to the image plane. The first row index is zero.
\delta_j - Row pixel resolution of the Pixel Spacing (0028,0030) attribute in units of mm.
Additional constraints apply:
1) The row and column direction cosine vectors shall be orthogonal, i.e. their dot product shall be zero.
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.
Example code for this is available here : http://www.dclunie.com/medical-image-faq/html/part2.html
Consider a DICOM file with direction cosines:
1.000000 0.000000 0.000000 0.000000 0.000000 -1.000000
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]
From : http://caddlab.rad.unc.edu/software/MetaIO/MetaIO-Introduction.pdf AnatomicalOrientation MET_STRING Specify anatomic ordering of axis. For example, RAS mean x-axis is Right to Left, y-axis is Anterior to Posterior, and Z-axis is Superior to Inferior. Use only [R|L] | [A|P] | [S|I] per axis.
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.
The basic concept of template recursion is the following
template <unsigned int N>
class Series
{
public:
void Compute()
{
sum += vector[N-1];
Series<N-1>::Compute(); // recursion
}
private:
int vector[N];
int sum;
}
void Series<1>::Compute() // Specialization for terminating the recursion
{
m_sum += vector[0];
}
The current proposal would involve to
- Move the TranformIndexToPhysicalPoint from itk::Image to itk::ImageBase
- Move the TranformIndexToContinuousIndex from itk::Image to itk::ImageBase
- Move the TranformPhysicalPointToIndex from itk::Image to itk::ImageBase
- Add an ImageHelper class (classname to be discussed)
- Invoke the Helper from the three methods above
A tentative version of the helper is in the following code:
template <unsigned int NImageDimension, class TPoint, unsigned int N>
class ImageHelper
{
public: typedef ImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::SpacingType SpacingType; typedef typename ImageType::PointType OriginType; typedef typename TPoint::ValueType TCoordRep;
itkStaticConstMacro(ImageDimension, unsigned int,
NImageDimension );
itkStaticConstMacro(DN, unsigned int, N-1);
static void TransformIndexToPhysicalPoint(
const SpacingType & spacing,
const OriginType & origin,
const IndexType & index, TPoint & point )
{
point[DN] = static_cast<TCoordRep>( spacing[DN] *
static_cast<double>( index[DN] ) + origin[DN] );
// Do next dimension, in order to follow the recursion
ImageHelper<NImageDimension,TPoint,DN>
::TransformIndexToPhysicalPoint(spacing,origin,index,point);
}
};
