Get voxel to RAS transformation from NIfTI file

TL;DR

I am struggling with representing points in different coordinate systems. I am using nilearn.

Is there a way to get the vox2ras_tkr and vox2ras transforms from a NIfTI file? Or can it only be obtained from an .mgz file?

Background

  • I have a set of points (anatomical landmarks: Nasion, left and right Preauricular Points [LPA/RPA])
  • I need to convert these points to be used with nilearn’s plot_anat function, so that I can plot the three points as shown in the code snippet below:
fig, axs = plt.subplots(3, 1)
for i in range(3):
    plot_anat(t1_fname, axes=axs[i], cut_coords=mri_pos[i, :])
plt.show()

However, to get the mri_pos points in the correct coordinate system, I need to run the following code:

# given `mri_pos`, obtained from visually aligning points in a GUI ...
t1 = nib.load(t1_fname)
vox2ras_tkr = t1.header.get_vox2ras_tkr()
ras2vox_tkr = np.linalg.inv(vox2ras_tkr)
vox2ras = t1.header.get_vox2ras()
mri_pos = apply_trans(ras2vox_tkr, mri_pos)  # in vox
mri_pos = apply_trans(vox2ras, mri_pos) # in RAS

Problem

  • When using a .mgz file as my t1_fname, everything works as expected
  • … when using a .nii.gz file however, I get the following error:

AttributeError: ‘Nifti1Header’ object has no attribute ‘get_vox2ras_tkr’

Is this to be expected?

Hi @sappelhoff, yes, not all image types have the same methods. While they share a common API, some have format-specific extensions, and get_vox2ras_tkr seems to be a specific transform of the normal affine matrix that’s relevant in FreeSurfer contexts.

Probably the simplest way to get this is:

import nibabel as nib
nii = nib.load(nii_fname)
mgh = nib.MGHImage(nii.dataobj, nii.affine)
vox2ras_tkr = mgh.header.get_vox2ras_tkr()

That said, is there an advantage to using get_vox2ras_tkr instead of working with img.affine directly? I would give that a shot. To the extent you can use the common API, your code can be less format-specific.

1 Like

Thanks @effigies

@agramfort has written another solution here. I paste it for convenience:

import numpy as np
import nibabel as nib

t1 = nib.load(<mypath>)


def get_vox2ras_tkr(t1):
    '''Get the vox2ras-tkr transform. Inspired
    by get_vox2ras_tkr in
    '''
    ds = t1.header._structarr['pixdim'][1:4]
    ns = t1.header._structarr['dim'][1:4] * ds / 2.0
    v2rtkr = np.array([[-ds[0], 0, 0, ns[0]],
                       [0, 0, ds[2], -ns[2]],
                       [0, -ds[1], 0, ns[1]],
                       [0, 0, 0, 1]], dtype=np.float32)
    return v2rtkr

np.testing.assert_array_equal(get_vox2ras_tkr(t1), t1_mgz.header.get_vox2ras_tkr())

it seems to me that img.affine corresponds to get_vox2ras … but not get_vox2ras_tkr? What does the tkr stand for?

tkr probably refers to tkregister.

And I agree that img.affine is not the same as get_vox2ras_tkr, but it may be that you can work with it and avoid all this mess. I don’t know what the advantage to working in the tkr variant is. It appears to be in LSI space with the centroid in the middle of the volume. If you’re not attached to that format, and your images are in the same space, you could probably just use img.affine.

The get_vox2ras_tkr function you’ve proposed assumes quite strongly (due to 'dim' and 'pixdim') that you’re working with Analyze-based images like NIfTI, which limits the applicability. The MGHImage approach I suggested will work with almost all nibabel-supported image types (I think GIFTI will fail and CIFTI-2 might sometimes fail to have a useful affine). If you want to avoid the indirection of building up an image object, you can also just pull the necessary info out of the zooms (deltas, in FreeSurfer terms) and shape, which are part of the SpatialImage API:

def get_vox2ras_tkr(t1):
    '''Get the vox2ras-tkr transform. Inspired
    by get_vox2ras_tkr in
    '''
    ds = t1.header.get_zooms()[:3]
    ns = np.array(t1.shape[:3]) * ds / 2.0
    v2rtkr = np.array([[-ds[0], 0, 0, ns[0]],
                       [0, 0, ds[2], -ns[2]],
                       [0, -ds[1], 0, ns[1]],
                       [0, 0, 0, 1]], dtype=np.float32)
    return v2rtkr