ijk2xyz: ijk2xyz
In AnalyzeFMRI: Functions for Analysis of fMRI Datasets Stored in the ANALYZE or NIFTI Format
Description
Usage
Arguments
Details
Value
See Also
Examples
Description
This function maps from data coordinates (e.g. column i, row j, slice k), into some real world (x,y,z) positions in space. These positions could relate to Talairach-Tournoux (T&T) space, MNI space, or patient-based scanner coordinates.
Usage
1
Arguments
ijk
matrix. Each column of ijk should contain a voxel index coordinates (i,j,k) to be mapped to its (x,y,z) real coordinates in some other space
method
1 (qform.code=sform.code=0),2 (qform.code>0, rigid
transformation) or 3 (sform.code>0, affine transformation).
L
header list of a NIFTI file
Details
The NIfTI format allows storage on disk to be in either a left- or
right-handed coordinate system. However, the format includes an implicit
spatial transformation into a RIGHT-HANDED coordinate system. This
transform maps from data coordinates (e.g. column i, row j,
slice k), into some real world (x,y,z) positions in
space. These positions could relate to Talairach-Tournoux (T&T) space,
MNI space, or patient-based scanner coordinates. For T&T, and MNI
coordinates, x increases from left to right, y increases
from posterior to anterior, and z increases in the inferior to
superior direction. Directions in the scanner coordinate system are
similar. MRI data is usually exported as DICOM format, which encodes the
positions and orientations of the slices. When data are converted from
DICOM to NIfTI-1 format, the relevant information can be determined from
the Pixel Spacing, Image Orientation (Patient) and Image Position
(Patient) fields of the DICOM files. NIfTI-1 also allows the space of one image to be mapped to that of
another (via a rigid or affine transform). This is to enable on-the-fly
resampling of registered images. This would allow intra-subject images,
collected with lots of different orientations or resolutions, to be
treated as if they are all in register.
Neurological and radiological conventions only relate to visualization
of axial images. They are unrelated to how the data are stored on disk,
or even how the real-world coordinates are represented.
It is more appropriate to consider whether the real-world coordinate
system is left- or right-handed (see below). Talairach and Tournoux use
a right-handed system, whereas the storage convention of ANALYZE files
is usually considered as left-handed. These coordinate systems are
mirror images of each other (if you are a psychologist, try explaining
why mirror images appear to be left-right flipped, rather than flipped
up-down, or back-front). Transforming between left- and right-handed
coordinate systems involves flipping, and can not be done by rotations
alone.
x=thumb, y=index finger (forefinger), z=left (resp. right) hand's middle finger for left-handed persons (resp. right-handed persons).
Volume orientation is given by a transformation that maps voxel indices
(i,j,k) to spatial coordinates (x,y,z),
typically anatomical coordinates assigned by the scanner. This
transformation (Method 2 in the ‘nifti1.h’ documentation)
is generated using the voxel dimensions, a quaternion encoding a rotation matrix, and a 3D shift, all stored in the NIfTI-1 header; details can be found in the ‘nifti1.h’ comments.
The NIfTI-1 header also provides for a general affine transformation, separate from that described by Method 2. This transformation (Method 3) also maps voxel indices (i,j,k) to (x,y,z), which in this case are typically coordinates in a standard space such as the Talairach space. The elements of this transformation matrix are stored in the NIfTI-1 header. For example, the Method 2 transformation can be constructed from the attributes from a set of DICOM files; the Method 3 transform can be computed offline and inserted into the header later.
The exact "meaning" of the coordinates given by the Method 2 and
Method 3 transformations is recorded in header fields
qform.code and sform.code, respectively. Code values can
indicate if the (x,y,z) axes are
Anatomical coordinates from the scanner (e.g., the DICOM header)
Aligned to some anatomical "truth" or standard
Aligned and warped to Talairach-Tournoux coordinates
Aligned and warped to MNI-152 coordinates
It is possible that neither transformation is specified (i.e.,
qform.code=sform.code=0), in which case we are left with the voxel size
in pixdim[], and no orientation is given or assumed. This use
(Method 1) is discouraged.
The basic idea behind having two coordinate systems is to allow the image to store information about (1) the scanner coordinate system used in the acquisition of the volume (in the qform) and (2) the relationship to a standard coordinate system - e.g. MNI coordinates (in the sform).
The qform allows orientation information to be kept for alignment purposes without losing volumetric information, since the qform only stores a rigid-body transformation which preserves volume. On the other hand, the sform stores a general affine transformation which can map the image coordinates into a standard coordinate system, like Talairach or MNI, without the need to resample the image.
By having both coordinate systems, it is possible to keep the original data (without resampling), along with information on how it was acquired (qform) and how it relates to other images via a standard space (sform). This ability is advantageous for many analysis pipelines, and has previously required storing additional files along with the image files. By using NIfTI-1 this extra information can be kept in the image files themselves.
Note: the qform and sform also store information on whether the coordinate system is left-handed or right-handed (see Q15) and so when both are set they must be consistent, otherwise the handedness of the coordinate system (often used to distinguish left-right order) is unknown and the results of applying operations to such an image are unspecified.
There are 3 different methods by which continuous coordinates can be
attached to voxels. The discussion below emphasizes 3D volumes, and
the continuous coordinates are referred to as (x,y,z). The voxel
index coordinates (i.e., the array indexes) are referred to as (i,j,k),
with valid ranges:
-
i = 0,…,dim[1]-1
-
j = 0,…,dim[2]-1 (if dim[0] >= 2)
-
k = 0,…,dim[3]-1 (if dim[0] >= 3)
The (x,y,z) coordinates refer to the CENTER of a voxel. In methods
2 and 3, the (x,y,z) axes refer to a subject-based coordinate system,
with
+x = Right +y = Anterior +z = Superior.
This is a right-handed coordinate system. However, the exact direction
these axes point with respect to the subject depends on qform.code
(Method 2) and sform.code (Method 3).
N.B.: The i index varies most rapidly, j index next, k index slowest.
Thus, voxel (i,j,k) is stored starting at location
(i + j*\code{dim[1]} + k*\code{dim[1]}*\code{dim[2]}) *
(\code{bitpix}/8)
into the dataset array.
N.B.: The ANALYZE 7.5 coordinate system is
+x = Left +y = Anterior +z = Superior
which is a left-handed coordinate system. This backwardness is
too difficult to tolerate, so this NIFTI-1 standard specifies the
coordinate order which is most common in functional neuroimaging.
N.B.: The 3 methods below all give the locations of the voxel centers
in the (x,y,z) coordinate system. In many cases, programs will wish
to display image data on some other grid. In such a case, the program
will need to convert its desired (x,y,z) values into (i,j,k) values
in order to extract (or interpolate) the image data. This operation
would be done with the inverse transformation to those described below.
N.B.: Method 2 uses a factor qfac which is either -1 or 1; qfac is
stored in the otherwise unused pixdim[0]. If pixdim[0]=0.0 (which
should not occur), we take qfac=1. Of course, pixdim[0] is only used
when reading a NIFTI-1 header, not when reading an ANALYZE 7.5 header.
N.B.: The units of (x,y,z) can be specified using the xyzt.units field.
-
METHOD 1 (the "old" way, used only when qform.code = 0):
The coordinate mapping from (i,j,k) to (x,y,z) is the ANALYZE
7.5 way. This is a simple scaling relationship:
x = pixdim[1] * i
y = pixdim[2] * j
z = pixdim[3] * k
No particular spatial orientation is attached to these (x,y,z)
coordinates. (NIFTI-1 does not have the ANALYZE 7.5 orient field,
which is not general and is often not set properly.) This method
is not recommended, and is present mainly for compatibility with
ANALYZE 7.5 files.
-
METHOD 2 (used when qform.code > 0, which should be the "normal" case):
The (x,y,z) coordinates are given by the pixdim[] scales, a rotation
matrix, and a shift. This method is intended to represent
"scanner-anatomical" coordinates, which are often embedded in the
image header (e.g., DICOM fields (0020,0032), (0020,0037), (0028,0030),
and (0018,0050)), and represent the nominal orientation and location of
the data. This method can also be used to represent "aligned"
coordinates, which would typically result from some post-acquisition
alignment of the volume to a standard orientation (e.g., the same
subject on another day, or a rigid rotation to true anatomical
orientation from the tilted position of the subject in the scanner).
The formula for (x,y,z) in terms of header parameters and (i,j,k) is:
[ x ] [ R11 R12 R13 ] [ pixdim[1] * i ] [ qoffset.x ]
[ y ] = [ R21 R22 R23 ] [ pixdim[2] * j ] + [ qoffset.y ]
[ z ] [ R31 R32 R33 ] [ qfac *
pixdim[3] * k ] [ qoffset.z ]
The qoffset.* shifts are in the NIFTI-1 header. Note that the center
of the (i,j,k)=(0,0,0) voxel (first value in the dataset array) is
just (x,y,z)=(\code{qoffset.x},\code{qoffset.y},\code{qoffset.z}).
The rotation matrix R is calculated from the quatern.* parameters.
This calculation is described below.
The scaling factor qfac is either 1 or -1. The rotation matrix R
defined by the quaternion parameters is "proper" (has determinant 1).
This may not fit the needs of the data; for example, if the image
grid is
i increases from Left-to-Right
j increases from Anterior-to-Posterior
k increases from Inferior-to-Superior
Then (i,j,k) is a left-handed triple. In this example, if qfac=1,
the R matrix would have to be
[ 1 0 0 ]
[ 0 -1 0 ] which is "improper" (determinant = -1).
[ 0 0 1 ]
If we set qfac=-1, then the R matrix would be
[ 1 0 0 ]
[ 0 -1 0 ] which is proper.
[ 0 0 -1 ]
This R matrix is represented by quaternion [a,b,c,d] = [0,1,0,0]
(which encodes a 180 degree rotation about the x-axis).
-
METHOD 3 (used when sform.code > 0):
The (x,y,z) coordinates are given by a general affine transformation
of the (i,j,k) indexes:
x = srow.x[0] * i + srow.x[1] * j + srow.x[2] * k + srow.x[3]
y = srow.y[0] * i + srow.y[1] * j + srow.y[2] * k + srow.y[3]
z = srow.z[0] * i + srow.z[1] * j + srow.z[2] * k + srow.z[3]
The srow.* vectors are in the NIFTI.1 header. Note that no use is
made of pixdim[] in this method.
WHY 3 METHODS?
Method 1 is provided only for backwards compatibility. The intention
is that Method 2 (qform.code > 0) represents the nominal voxel locations
as reported by the scanner, or as rotated to some fiducial orientation and
location. Method 3, if present (sform.code > 0), is to be used to give
the location of the voxels in some standard space. The sform.code
indicates which standard space is present. Both methods 2 and 3 can be
present, and be useful in different contexts (method 2 for displaying the
data on its original grid; method 3 for displaying it on a standard grid).
In this scheme, a dataset would originally be set up so that the
Method 2 coordinates represent what the scanner reported. Later,
a registration to some standard space can be computed and inserted
in the header. Image display software can use either transform,
depending on its purposes and needs.
In Method 2, the origin of coordinates would generally be whatever
the scanner origin is; for example, in MRI, (0,0,0) is the center
of the gradient coil.
In Method 3, the origin of coordinates would depend on the value
of sform.code; for example, for the Talairach coordinate system,
(0,0,0) corresponds to the Anterior Commissure.
QUATERNION REPRESENTATION OF ROTATION MATRIX (METHOD 2)
The orientation of the (x,y,z) axes relative to the (i,j,k) axes
in 3D space is specified using a unit quaternion [a,b,c,d], where
a*a+b*b+c*c+d*d=1. The (b,c,d) values are all that is needed, since
we require that a = √(1.0-(b*b+c*c+d*d)) be nonnegative. The (b,c,d)
values are stored in the (quatern.b,quatern.c,quatern.d) fields.
The quaternion representation is chosen for its compactness in
representing rotations. The (proper) 3x3 rotation matrix that
corresponds to [a,b,c,d] is
[ a*a+b*b-c*c-d*d 2*b*c-2*a*d 2*b*d+2*a*c ]
R = [ 2*b*c+2*a*d a*a+c*c-b*b-d*d 2*c*d-2*a*b ]
[ 2*b*d-2*a*c 2*c*d+2*a*b a*a+d*d-c*c-b*b ]
[ R11 R12 R13 ]
= [ R21 R22 R23 ]
[ R31 R32 R33 ]
If (p,q,r) is a unit 3-vector, then rotation of angle h about that
direction is represented by the quaternion
[a,b,c,d] = [cos(h/2), p*sin(h/2), q*sin(h/2), r*sin(h/2)].
Requiring a >= 0 is equivalent to requiring -Pi <= h <= Pi. (Note that
[-a,-b,-c,-d] represents the same rotation as [a,b,c,d]; there are 2
quaternions that can be used to represent a given rotation matrix R.)
To rotate a 3-vector (x,y,z) using quaternions, we compute the
quaternion product
[0,x',y',z'] = [a,b,c,d] * [0,x,y,z] * [a,-b,-c,-d]
which is equivalent to the matrix-vector multiply
[ x' ] [ x ]
[ y' ] = R [ y ] (equivalence depends on a*a+b*b+c*c+d*d=1)
[ z' ] [ z ]
Multiplication of 2 quaternions is defined by the following:
[a,b,c,d] = a*1 + b*I + c*J + d*K
where
I*I = J*J = K*K = -1 (I,J,K are square roots of -1)
I*J = K , J*K = I , K*I = J
J*I = -K , K*J = -I , I*K = -J (not
commutative!).
For example
[a,b,0,0] * [0,0,0,1] = [0,0,-b,a]
since this expands to
(a+b*I)*(K) = (a*K+b*I*K) = (a*K-b*J).
The above formula shows how to go from quaternion (b,c,d) to
rotation matrix and direction cosines. Conversely, given R,
we can compute the fields for the NIFTI-1 header by
a = 0.5 * √(1+R11+R22+R33) (not stored)
b = 0.25 * (R32-R23) / a => quatern.b
c = 0.25 * (R13-R31) / a => quatern.c
d = 0.25 * (R21-R12) / a => quatern.d
If a=0 (a 180 degree rotation), alternative formulas are needed.
See the ‘nifti1.io.c’ function mat44.to.quatern() for an implementation
of the various cases in converting R to [a,b,c,d].
Note that R-transpose (= R-inverse) would lead to the quaternion
[a,-b,-c,-d].
The choice to specify the qoffset.x (etc.) values in the final
coordinate system is partly to make it easy to convert DICOM images to
this format. The DICOM attribute "Image Position (Patient)" (0020,0032)
stores the (Xd,Yd,Zd) coordinates of the center of the first voxel.
Here, (Xd,Yd,Zd) refer to DICOM coordinates, and Xd=-x, Yd=-y, Zd=z,
where (x,y,z) refers to the NIFTI coordinate system discussed above.
(i.e., DICOM +Xd is Left, +Yd is Posterior, +Zd is Superior,
whereas +x is Right, +y is Anterior , +z is Superior. )
Thus, if the (0020,0032) DICOM attribute is extracted into (px,py,pz), then
qoffset.x = -px qoffset.y = -py qoffset.z = pz
is a reasonable setting when qform.code=NIFTI.XFORM.SCANNER.ANAT.
That is, DICOM's coordinate system is 180 degrees rotated about the z-axis
from the neuroscience/NIFTI coordinate system. To transform between DICOM
and NIFTI, you just have to negate the x- and y-coordinates.
The DICOM attribute (0020,0037) "Image Orientation (Patient)" gives the
orientation of the x- and y-axes of the image data in terms of 2 3-vectors.
The first vector is a unit vector along the x-axis, and the second is
along the y-axis. If the (0020,0037) attribute is extracted into the
value (xa,xb,xc,ya,yb,yc), then the first two columns of the R matrix
would be
[ -xa -ya ]
[ -xb -yb ]
[ xc yc ]
The negations are because DICOM's x- and y-axes are reversed relative
to NIFTI's. The third column of the R matrix gives the direction of
displacement (relative to the subject) along the slice-wise direction.
This orientation is not encoded in the DICOM standard in a simple way;
DICOM is mostly concerned with 2D images. The third column of R will be
either the cross-product of the first 2 columns or its negative. It is
possible to infer the sign of the 3rd column by examining the coordinates
in DICOM attribute (0020,0032) "Image Position (Patient)" for successive
slices. However, this method occasionally fails for reasons that I
(RW Cox) do not understand.
Value
A list containing the matrix xyz of the positions of the points
specified in ijk.
See Also
Examples
1
2
3
4
L <- f.read.header(system.file("example-nifti.hdr",
package="AnalyzeFMRI"))
ijk <- matrix(c(1,1,1,2,3,7),byrow=FALSE,nrow=3)
ijk2xyz(ijk=ijk,method=2,L)
Example output
Loading required package: R.matlab
R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.
Attaching package: 'R.matlab'
The following objects are masked from 'package:base':
getOption, isOpen
Loading required package: fastICA
Loading required package: tcltk
Loading required package: tkrplot
Warning messages:
1: no DISPLAY variable so Tk is not available
2: loading Rplot failed
$xyz
[,1] [,2]
[1,] 78 75
[2,] -111 -105
[3,] -51 -33
AnalyzeFMRI documentation built on Oct. 5, 2021, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value See Also Examples
Description
This function maps from data coordinates (e.g. column i, row j, slice k), into some real world (x,y,z) positions in space. These positions could relate to Talairach-Tournoux (T&T) space, MNI space, or patient-based scanner coordinates.
Usage
1 |
Arguments
ijk |
matrix. Each column of ijk should contain a voxel index coordinates (i,j,k) to be mapped to its (x,y,z) real coordinates in some other space |
method |
1 (qform.code=sform.code=0),2 (qform.code>0, rigid transformation) or 3 (sform.code>0, affine transformation). |
L |
header list of a NIFTI file |
Details
The NIfTI format allows storage on disk to be in either a left- or right-handed coordinate system. However, the format includes an implicit spatial transformation into a RIGHT-HANDED coordinate system. This transform maps from data coordinates (e.g. column i, row j, slice k), into some real world (x,y,z) positions in space. These positions could relate to Talairach-Tournoux (T&T) space, MNI space, or patient-based scanner coordinates. For T&T, and MNI coordinates, x increases from left to right, y increases from posterior to anterior, and z increases in the inferior to superior direction. Directions in the scanner coordinate system are similar. MRI data is usually exported as DICOM format, which encodes the positions and orientations of the slices. When data are converted from DICOM to NIfTI-1 format, the relevant information can be determined from the Pixel Spacing, Image Orientation (Patient) and Image Position (Patient) fields of the DICOM files. NIfTI-1 also allows the space of one image to be mapped to that of another (via a rigid or affine transform). This is to enable on-the-fly resampling of registered images. This would allow intra-subject images, collected with lots of different orientations or resolutions, to be treated as if they are all in register.
Neurological and radiological conventions only relate to visualization of axial images. They are unrelated to how the data are stored on disk, or even how the real-world coordinates are represented. It is more appropriate to consider whether the real-world coordinate system is left- or right-handed (see below). Talairach and Tournoux use a right-handed system, whereas the storage convention of ANALYZE files is usually considered as left-handed. These coordinate systems are mirror images of each other (if you are a psychologist, try explaining why mirror images appear to be left-right flipped, rather than flipped up-down, or back-front). Transforming between left- and right-handed coordinate systems involves flipping, and can not be done by rotations alone.
x=thumb, y=index finger (forefinger), z=left (resp. right) hand's middle finger for left-handed persons (resp. right-handed persons).
Volume orientation is given by a transformation that maps voxel indices
(i,j,k) to spatial coordinates (x,y,z),
typically anatomical coordinates assigned by the scanner. This
transformation (Method 2 in the ‘nifti1.h’ documentation)
is generated using the voxel dimensions, a quaternion encoding a rotation matrix, and a 3D shift, all stored in the NIfTI-1 header; details can be found in the ‘nifti1.h’ comments.
The NIfTI-1 header also provides for a general affine transformation, separate from that described by Method 2. This transformation (Method 3) also maps voxel indices (i,j,k) to (x,y,z), which in this case are typically coordinates in a standard space such as the Talairach space. The elements of this transformation matrix are stored in the NIfTI-1 header. For example, the Method 2 transformation can be constructed from the attributes from a set of DICOM files; the Method 3 transform can be computed offline and inserted into the header later.
The exact "meaning" of the coordinates given by the Method 2 and
Method 3 transformations is recorded in header fields
qform.code and sform.code, respectively. Code values can
indicate if the (x,y,z) axes are
Anatomical coordinates from the scanner (e.g., the DICOM header)
Aligned to some anatomical "truth" or standard
Aligned and warped to Talairach-Tournoux coordinates
Aligned and warped to MNI-152 coordinates
It is possible that neither transformation is specified (i.e.,
qform.code=sform.code=0), in which case we are left with the voxel size
in pixdim[], and no orientation is given or assumed. This use
(Method 1) is discouraged.
The basic idea behind having two coordinate systems is to allow the image to store information about (1) the scanner coordinate system used in the acquisition of the volume (in the qform) and (2) the relationship to a standard coordinate system - e.g. MNI coordinates (in the sform).
The qform allows orientation information to be kept for alignment purposes without losing volumetric information, since the qform only stores a rigid-body transformation which preserves volume. On the other hand, the sform stores a general affine transformation which can map the image coordinates into a standard coordinate system, like Talairach or MNI, without the need to resample the image.
By having both coordinate systems, it is possible to keep the original data (without resampling), along with information on how it was acquired (qform) and how it relates to other images via a standard space (sform). This ability is advantageous for many analysis pipelines, and has previously required storing additional files along with the image files. By using NIfTI-1 this extra information can be kept in the image files themselves.
Note: the qform and sform also store information on whether the coordinate system is left-handed or right-handed (see Q15) and so when both are set they must be consistent, otherwise the handedness of the coordinate system (often used to distinguish left-right order) is unknown and the results of applying operations to such an image are unspecified.
There are 3 different methods by which continuous coordinates can be
attached to voxels. The discussion below emphasizes 3D volumes, and
the continuous coordinates are referred to as (x,y,z). The voxel
index coordinates (i.e., the array indexes) are referred to as (i,j,k),
with valid ranges:
-
i = 0,…,
dim[1]-1
-
j = 0,…,
dim[2]-1(ifdim[0]>= 2)
-
k = 0,…,
dim[3]-1(ifdim[0]>= 3)
The (x,y,z) coordinates refer to the CENTER of a voxel. In methods
2 and 3, the (x,y,z) axes refer to a subject-based coordinate system,
with
+x = Right +y = Anterior +z = Superior.
This is a right-handed coordinate system. However, the exact direction
these axes point with respect to the subject depends on qform.code
(Method 2) and sform.code (Method 3).
N.B.: The i index varies most rapidly, j index next, k index slowest. Thus, voxel (i,j,k) is stored starting at location
(i + j*\code{dim[1]} + k*\code{dim[1]}*\code{dim[2]}) * (\code{bitpix}/8)
into the dataset array.
N.B.: The ANALYZE 7.5 coordinate system is
+x = Left +y = Anterior +z = Superior
which is a left-handed coordinate system. This backwardness is too difficult to tolerate, so this NIFTI-1 standard specifies the coordinate order which is most common in functional neuroimaging.
N.B.: The 3 methods below all give the locations of the voxel centers in the (x,y,z) coordinate system. In many cases, programs will wish to display image data on some other grid. In such a case, the program will need to convert its desired (x,y,z) values into (i,j,k) values in order to extract (or interpolate) the image data. This operation would be done with the inverse transformation to those described below.
N.B.: Method 2 uses a factor qfac which is either -1 or 1; qfac is
stored in the otherwise unused pixdim[0]. If pixdim[0]=0.0 (which
should not occur), we take qfac=1. Of course, pixdim[0] is only used
when reading a NIFTI-1 header, not when reading an ANALYZE 7.5 header.
N.B.: The units of (x,y,z) can be specified using the xyzt.units field.
-
METHOD 1 (the "old" way, used only when
qform.code= 0):The coordinate mapping from (i,j,k) to (x,y,z) is the ANALYZE 7.5 way. This is a simple scaling relationship:
x = pixdim[1] * i
y = pixdim[2] * j
z = pixdim[3] * kNo particular spatial orientation is attached to these (x,y,z) coordinates. (NIFTI-1 does not have the ANALYZE 7.5 orient field, which is not general and is often not set properly.) This method is not recommended, and is present mainly for compatibility with ANALYZE 7.5 files.
-
METHOD 2 (used when
qform.code> 0, which should be the "normal" case):The (x,y,z) coordinates are given by the
pixdim[]scales, a rotation matrix, and a shift. This method is intended to represent "scanner-anatomical" coordinates, which are often embedded in the image header (e.g., DICOM fields (0020,0032), (0020,0037), (0028,0030), and (0018,0050)), and represent the nominal orientation and location of the data. This method can also be used to represent "aligned" coordinates, which would typically result from some post-acquisition alignment of the volume to a standard orientation (e.g., the same subject on another day, or a rigid rotation to true anatomical orientation from the tilted position of the subject in the scanner). The formula for (x,y,z) in terms of header parameters and (i,j,k) is:[ x ] [ R11 R12 R13 ] [ pixdim[1]* i ][ qoffset.x][ y ] = [ R21 R22 R23 ] [ pixdim[2]* j ]+ [ qoffset.y][ z ] [ R31 R32 R33 ] [ qfac*pixdim[3]* k ][ qoffset.z]The
qoffset.* shifts are in the NIFTI-1 header. Note that the center of the (i,j,k)=(0,0,0) voxel (first value in the dataset array) is just (x,y,z)=(\code{qoffset.x},\code{qoffset.y},\code{qoffset.z}).The rotation matrix R is calculated from the
quatern.* parameters. This calculation is described below.The scaling factor
qfacis either 1 or -1. The rotation matrix R defined by the quaternion parameters is "proper" (has determinant 1). This may not fit the needs of the data; for example, if the image grid is
i increases from Left-to-Right
j increases from Anterior-to-Posterior
k increases from Inferior-to-Superior
Then (i,j,k) is a left-handed triple. In this example, ifqfac=1, the R matrix would have to be[ 1 0 0 ]
[ 0 -1 0 ] which is "improper" (determinant = -1).
[ 0 0 1 ]If we set
qfac=-1, then the R matrix would be[ 1 0 0 ]
[ 0 -1 0 ] which is proper.
[ 0 0 -1 ]This R matrix is represented by quaternion [a,b,c,d] = [0,1,0,0] (which encodes a 180 degree rotation about the x-axis).
-
METHOD 3 (used when
sform.code> 0):The (x,y,z) coordinates are given by a general affine transformation of the (i,j,k) indexes:
x = srow.x[0]* i +srow.x[1]* j +srow.x[2]* k +srow.x[3]y = srow.y[0]* i +srow.y[1]* j +srow.y[2]* k +srow.y[3]z = srow.z[0]* i +srow.z[1]* j +srow.z[2]* k +srow.z[3]The
srow.* vectors are in the NIFTI.1 header. Note that no use is made ofpixdim[]in this method. WHY 3 METHODS?
Method 1 is provided only for backwards compatibility. The intention is that Method 2 (
qform.code> 0) represents the nominal voxel locations as reported by the scanner, or as rotated to some fiducial orientation and location. Method 3, if present (sform.code> 0), is to be used to give the location of the voxels in some standard space. Thesform.codeindicates which standard space is present. Both methods 2 and 3 can be present, and be useful in different contexts (method 2 for displaying the data on its original grid; method 3 for displaying it on a standard grid).In this scheme, a dataset would originally be set up so that the Method 2 coordinates represent what the scanner reported. Later, a registration to some standard space can be computed and inserted in the header. Image display software can use either transform, depending on its purposes and needs.
In Method 2, the origin of coordinates would generally be whatever the scanner origin is; for example, in MRI, (0,0,0) is the center of the gradient coil.
In Method 3, the origin of coordinates would depend on the value of
sform.code; for example, for the Talairach coordinate system, (0,0,0) corresponds to the Anterior Commissure.QUATERNION REPRESENTATION OF ROTATION MATRIX (METHOD 2)
The orientation of the (x,y,z) axes relative to the (i,j,k) axes in 3D space is specified using a unit quaternion [a,b,c,d], where a*a+b*b+c*c+d*d=1. The (b,c,d) values are all that is needed, since we require that a = √(1.0-(b*b+c*c+d*d)) be nonnegative. The (b,c,d) values are stored in the (
quatern.b,quatern.c,quatern.d) fields.The quaternion representation is chosen for its compactness in representing rotations. The (proper) 3x3 rotation matrix that corresponds to [a,b,c,d] is
[ a*a+b*b-c*c-d*d 2*b*c-2*a*d 2*b*d+2*a*c ] R = [ 2*b*c+2*a*d a*a+c*c-b*b-d*d 2*c*d-2*a*b ] [ 2*b*d-2*a*c 2*c*d+2*a*b a*a+d*d-c*c-b*b ] [ R11 R12 R13 ] = [ R21 R22 R23 ] [ R31 R32 R33 ] If (p,q,r) is a unit 3-vector, then rotation of angle h about that direction is represented by the quaternion
[a,b,c,d] = [cos(h/2), p*sin(h/2), q*sin(h/2), r*sin(h/2)].
Requiring a >= 0 is equivalent to requiring -Pi <= h <= Pi. (Note that [-a,-b,-c,-d] represents the same rotation as [a,b,c,d]; there are 2 quaternions that can be used to represent a given rotation matrix R.) To rotate a 3-vector (x,y,z) using quaternions, we compute the quaternion product
[0,x',y',z'] = [a,b,c,d] * [0,x,y,z] * [a,-b,-c,-d]
which is equivalent to the matrix-vector multiply
[ x' ] [ x ] [ y' ] = R [ y ] (equivalence depends on a*a+b*b+c*c+d*d=1) [ z' ] [ z ] Multiplication of 2 quaternions is defined by the following:
[a,b,c,d] = a*1 + b*I + c*J + d*K
where
I*I = J*J = K*K = -1 (I,J,K are square roots of -1)
I*J = K , J*K = I , K*I = J
J*I = -K , K*J = -I , I*K = -J (not commutative!).For example
[a,b,0,0] * [0,0,0,1] = [0,0,-b,a]
since this expands to
(a+b*I)*(K) = (a*K+b*I*K) = (a*K-b*J).
The above formula shows how to go from quaternion (b,c,d) to rotation matrix and direction cosines. Conversely, given R, we can compute the fields for the NIFTI-1 header by
a = 0.5 * √(1+R11+R22+R33) (not stored)
b = 0.25 * (R32-R23) / a =>quatern.b
c = 0.25 * (R13-R31) / a =>quatern.c
d = 0.25 * (R21-R12) / a =>quatern.dIf
a=0(a 180 degree rotation), alternative formulas are needed. See the ‘nifti1.io.c’ function mat44.to.quatern() for an implementation of the various cases in converting R to [a,b,c,d].Note that R-transpose (= R-inverse) would lead to the quaternion [a,-b,-c,-d].
The choice to specify the
qoffset.x(etc.) values in the final coordinate system is partly to make it easy to convert DICOM images to this format. The DICOM attribute "Image Position (Patient)" (0020,0032) stores the (Xd,Yd,Zd) coordinates of the center of the first voxel. Here, (Xd,Yd,Zd) refer to DICOM coordinates, and Xd=-x, Yd=-y, Zd=z, where (x,y,z) refers to the NIFTI coordinate system discussed above. (i.e., DICOM +Xd is Left, +Yd is Posterior, +Zd is Superior,
whereas +x is Right, +y is Anterior , +z is Superior. )
Thus, if the (0020,0032) DICOM attribute is extracted into (px,py,pz), then
qoffset.x= -pxqoffset.y= -pyqoffset.z= pz
is a reasonable setting whenqform.code=NIFTI.XFORM.SCANNER.ANAT.That is, DICOM's coordinate system is 180 degrees rotated about the z-axis from the neuroscience/NIFTI coordinate system. To transform between DICOM and NIFTI, you just have to negate the x- and y-coordinates.
The DICOM attribute (0020,0037) "Image Orientation (Patient)" gives the orientation of the x- and y-axes of the image data in terms of 2 3-vectors. The first vector is a unit vector along the x-axis, and the second is along the y-axis. If the (0020,0037) attribute is extracted into the value (xa,xb,xc,ya,yb,yc), then the first two columns of the R matrix would be
[ -xa -ya ]
[ -xb -yb ]
[ xc yc ]
The negations are because DICOM's x- and y-axes are reversed relative to NIFTI's. The third column of the R matrix gives the direction of displacement (relative to the subject) along the slice-wise direction. This orientation is not encoded in the DICOM standard in a simple way; DICOM is mostly concerned with 2D images. The third column of R will be either the cross-product of the first 2 columns or its negative. It is possible to infer the sign of the 3rd column by examining the coordinates in DICOM attribute (0020,0032) "Image Position (Patient)" for successive slices. However, this method occasionally fails for reasons that I (RW Cox) do not understand.
Value
A list containing the matrix xyz of the positions of the points specified in ijk.
See Also
Examples
1 2 3 4 | L <- f.read.header(system.file("example-nifti.hdr",
package="AnalyzeFMRI"))
ijk <- matrix(c(1,1,1,2,3,7),byrow=FALSE,nrow=3)
ijk2xyz(ijk=ijk,method=2,L)
|
Example output
Loading required package: R.matlab
R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.
Attaching package: 'R.matlab'
The following objects are masked from 'package:base':
getOption, isOpen
Loading required package: fastICA
Loading required package: tcltk
Loading required package: tkrplot
Warning messages:
1: no DISPLAY variable so Tk is not available
2: loading Rplot failed
$xyz
[,1] [,2]
[1,] 78 75
[2,] -111 -105
[3,] -51 -33
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.