A colleague uses FSL FEAT to process fMRI data. I would like to use transformation matrix files produced by FLIRT during this process, such as standard2highres.mat, in 3D Slicer, but I do not understand how .mat affine matrices correspond to Slicer transformation matrices. In Slicer, affine transformation matrices are used to describe registration transformations between 3D physical spaces in RAS (right-anterior-superior) coordinate systems. FLIRT .mat files, I think, are similar, but it is possible that they describe transformations between IJK (image voxel coordinate) spaces. In any case, a FLIRT-produced .mat affine matrix and a Slicer-produced affine matrix are not the same for the same pair of images and the same quality of registration. I understand the Slicer system quite well, so I don’t have questions there, but I don’t understand the FLIRT conventions well enough to translate a FLIRT affine matrix to the equivalent Slicer version (which is what I want to use).
Could anyone here give me or point me to a description of what exactly the FLIRT affine matrix is doing and how it would be used? For example, if multiplied by a vector to produce a transformed point, what coordinate systems are the input and output points in? If v_out = F * v_in, where F is the FLIRT affine transformation matrix, is v_in:
a vector with coordinates in the moving image index space (e.g. the (i,j,k) voxel of the moving image array),
a vector with coordinate in the moving image physical space (e.g. for voxel (i,j,k) and voxel size ri x rj x rk (each in mm), the spatial coordinate would be (i * ri, j * rj, k * rk))
a vector with coordinate in the moving image physical space but specified as a location relative to the moving image center of mass (I’m not sure exactly how this would work, but I see notes that FLIRT transformation components are constructed such that the center of rotation and scaling is the image center of mass, though I don’t think this should actually matter for the final transformation matrix (but maybe I’m wrong))
something else?
I would appreciate if anyone here is able to shed any light on this subject. Thanks so much!
I think the FSL file is a simple voxel-to-voxel 3d transformation matrix. However, you need to be aware that some (but not all) FSL tools will mirror the column data (NIfTI i dimension) when loading a NIfTI image with a positive determinant (such that the loaded image always has a negative determinant). See my notes on the different behavior of fslstats and tfceS. So you will want to test your solution with source and target images with all permutations of positive and negative determinants.
Thanks so much for the reply, @neurolabusc, I appreciate it and will take your suggestion. Thanks also for the warning about different behavior depending on the determinant, I will watch out for that. I think I may be close to figuring this out, and will post back here if I think I am successful.
It is possible to convert between Slicer registration transform matrices and FSL/FLIRT registration transform matrices. They differ because Slicer transforms describe transformations within RAS space, whereas FLIRT transforms describe transformations between scaled (and sometimes reflected) voxel spaces. Converting between them requires understanding exactly when FLIRT transforms performed a reflection/translation before registration, as well as what other elements of the images are included in the FLIRT transformation matrix which are not part of the Slicer transform matrix. This is not well documented and was not intuitive to me (hence the long explanation below). If you want to skip to something which shows some details in code, take a look at the gist linked at the bottom of the post. That code is Slicer-based because that’s the tool I know well, but might be useful even if you were coming from the FSL side.
Detailed Explanation
Slicer transformation matrices relate RAS spatial coordinates. That is, image data for each image is first situated in the RAS spatial coordinate system, and then the Slicer transformation matrix converts a vector of RAS coordinate (for example the location of some anatomical feature in the moving image) to the corresponding RAS coordinate relative to the fixed image (the location of the same anatomical feature in the fixed image in spatial RAS space)
FSL/FLIRT transformation matrices relate scaled (and sometimes reflected) ijk spatial coordinates. The units are millimeters, but the header information regarding the orientation and location of image volume axes in a world coordinate system (the qform or sform in NIFTI files, or the IJKToRASMatrix in Slicer) is almost entirely ignored. I say “almost entirely” because FSL/FLIRT does look at the IJKToRASMatrix (or equivalent), and if the determinant of that matrix is positive, before registration is carried out, the orientation of the i-axis is reversed and the origin is moved to the other corner voxel on the i-axis. This reflection-translation inverts the sign of the determinant so that it is now negative.
It is possible to obtain the Slicer transformation matrix from the FSL/FLIRT transformation matrix, or vice versa, as long as you also have the two image volumes which were used in the registration which generated the matrix you have, or at the absolute minimum, image volumes which you know to have identical header data (same IJKToRASDirectionMatrix, voxel size, image array dimensions; the voxel contents are not needed) to the image volumes which were registered.
This was not straightforward to figure out, but I think I have it all worked out now, and am going to try to lay it out clearly enough that hopefully this can be helpful to others. In the final understanding, I found the series of coordinate systems in the table below very helpful. Some nomenclature notes first:
I use V for a generic image volume, F for the image volume used as the fixed image in the registration pair (“-ref” in the FLIRT command), and M for the moving image volume in the registration pair (“-in” in the FLIRT command). I use ijk to refer to image voxel array coordinates. I use hjk to refer to the image voxel array coordinates after the i-axis has been optionally flipped and shifted. Note that hjk is sometimes identical to ijk (if the determinant of the IJKToRasDirectionMatrix is negative), that is, I use hjk notation after the flip could have been applied, even if it wasn’t actually applied (because the determinant was already negative). I use hjkMm to refer the coordinate system where hjk indices have been scaled by the voxel size and therefore the coordinate system has millimeter units rather than voxel units. The FLIRT transformation matrix transforms between the moving image’s hjkMm coordinate system to the fixed image’s hjkMm coordinate system. The next coordinate system has the HJKToRASDirectionMatrix rotation applied to hjkMm coordinates; I call this one hras0, where the “h” is a reminder that the voxel origin may have been switched, the “ras” is because we are in a RAS coordinate system, and the “0” is a reminder that the ”h” voxel origin has not yet been moved from [0,0,0]. That is, the transformations from hjk to hras0 involve only scaling and rotation around the origin, so [0,0,0] in hjk space is also [0,0,0] in hras0 space. The last coordinate system is called ras and represents the completion of the transformation from ijk voxel array space to a spatial location in RAS space as defined by the image header. ras differs from hras0 only by a translation transformation which moves the location of [0,0,0] in ijk space, as represented in hras0 space, to the origin location specified in the image header. This ras coordinate system is the same one used as the 3D space RAS coordinate system in 3D Slicer. It no longer has an “h” in the notation because no element remains tying it to the h-axis flip/translation. Between ijk and ras, if a flip/translation occurred, it happened between the ijk coordinate system and the hjk coordinate system, and then was undone by the combination of moving from the hjkMm space to the hras0 space and the translation from hras0 space to ras.
Coord System Name
Units
Description
V_ijk
voxels
Voxel array coordinate system for volume V
V_hjk
voxels
Voxel array coordinate system for volume V, except that a h coordinate replaces the i coord, and where, Option A: h = i if the determinant of the IJKToRAS matrix for V is negative, or Option B: h = N - 1 - i if the determinant is negative. N = the number of image array voxels in the i dimension. This effectively switches the origin of the image from one corner of the array to the opposite corner (from the original [0,0,0] voxel to the original [N-1,0,0] voxel). This is done by FSL to enforce that all images have a fixed handedness (i.e. forcing a negative determinant). I’m not totally sure why anyone would do this; perhaps because if they didn’t they would have to allow reflections in their transform matrices?
V_hjkMm
mm
scaled voxel coordinate system. Origin has not moved relative to the V_hjk coordinate system, but [1,1,1] in the V_hjk coordinate system is located at [Sp_i, Sp_j, Sp_k] in the V_hjkMm coordinate system, where Sp_i is the voxel spacing in the i direction (note Sp_h and Sp_i are the same), and likewise for Sp_j, and Sp_k. More generally, the point [h,j,k] is located at [h*Sp_i, j*Sp_j, k*Sp_k] in the V_hjkMm coordinate system.
V_hras0
mm
rotated coordinate system. Origin has not moved relative to the V_hjk coordinate system, but the axes which point in the h, j, and k directions in the V_hjk and V_hjkMm coordinate systems now point in the directions of the columns of the V_HJKDirs matrix. The V_HJKDirs matrix is the same as the IJKToRAS direction matrix of V if det(V_IJKToRASDirs) < 0, but if det(V_IJKToRASDirs) > 0, then V_HJKDirs = V_IJKToRASDirs @ diag([-1, 1, 1, 1]), because the h axis is reversed with respect to the i axis. Note that V_IJKToRASDirs is IJKToRAS **DIRECTIONS** matrix, not the IJKToRASMatrix (which includes scaling and origin translation). In contrast, the IJKToRASDirection matrix (obtained by volNode.GetIJKToRASDirectionMatrix() in Slicer) has unit length columns and no translation component.
V_ras
mm
spatial coordinates in the Slicer world coordinate system. The axis directions match those in V_hras0, but there is a translation-only difference between the two coordinate systems. In V_hras0, the V_hjk origin voxel is still at [0,0,0], but in V_ras, that voxel is located at V_hOri, and the V_ijk origin voxel is located at volNode.GetOrigin(). V_hOri is calculated to compensate for the switch of origin from i to h (if it happened; if det(V_IJKToRASDirs)<0, h=i, and hOri=ori).
Here is a table of what each of the stepwise transformation matrices are between adjacent coordinate systems:
Transformation
Calculation
V_ijk_to_hjk
Identity matrix if det(IJKToRas) < 0, or translationMatrix([N-1, 0 , 0]) @ diag([-1, 1,1,1]) if det(IJKToRas) > 0, where @ means matrix multiplication (as in numpy)
V_hjk_to_hjkMm
diag([Sp_i, Sp_j, Sp_k, 1]) if the voxel spacing is [Sp_i, Sp_j, Sp_k] in the i, j, and k directions, respectively
V_hjkMm_to_hras0
IJKToRASDirectionMatrix if det(IJKToRas) < 0, or IJKToRASDirectionMatrix @ diag([-1, 1, 1, 1]) if det(IJKToRas) > 0
V_hras0_to_ras
translationMatrix(ori - ijk_origin_in_hras0), where ori is the origin listed in the header information for image volume V, and ijk_origin_in_hras0 is the location of the point with ijk=[0,0,0] in the hras0 coordinate system. This can be calculated as V_hjkMm_to_hras0 @ V_hjk_to_hjkMm @ V_ijk_to_hjk @ [0,0,0,1].
If a registration has been carried out between a moving image volume M and a fixed image volume F, then one can find the corresponding ijk coordinate in F to an ijk coordinate in M through the following chain of coordinate system transformations:
M_ijk to M_hjk to M_hjkMm to M_hras0 to M_ras to F_ras to F_hras0 to F_hjkMm to F_hjk to F_ijk
Slicer and FLIRT registrations respresent different portions of this chain. A Slicer registration transfrom matrix would represent just the M_ras to F_ras transformation, whereas a FLIRT registration transform matrix represents the M_hjkMm to M_hras0 to M_ras to F_ras to F_hras0 to F_hjkMm transformation.
With this knowledge, we can convert between them!
FLIRT transformation matrix from Slicer transformation matrix: