FreeSurfer coordinate spaces and translation

[effigies]

Following are some notes on handling "volume geometry" in FreeSurfer, and calculating the tkregister affine, the scanner affine, and the transform between the two.

This was motivated by a desire to understand what mris_convert --to-scanner is actually doing. Various pieces of this exist in nibabel, nitransforms, HCP pipelines and sMRIPrep, but there isn't really one place that puts it all together. Probably the closest is https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems, but this brings the underlying data structure to the fore and hopefully this will be a useful reference for those more comfortable reading Python.

Volume geometry

FreeSurfer encodes the affine as zooms, direction cosines, and c_ras, which is the the RAS coordinate of the center of the data array. We will model this with a simple dataclass:

from dataclass import dataclass

@dataclass
class VolGeom:
    # Using tuples to show shape. Will mostly be numpy arrays.
    shape: tuple[int, int, int]
    zooms: tuple[float, float, float]
    cosines: tuple[
        tuple[float, float, float],
        tuple[float, float, float],
        tuple[float, float, float]
    ]
    c_ras: tuple[float, float, float]

This volume geometry gets represented in an MGHHeader as follows:

>>> print(mgh.header)
<class 'nibabel.freesurfer.mghformat.MGHHeader'> object, endian='>'
version      : 1
dims         : [320 320 320   1]
type         : 0
dof          : 0
goodRASFlag  : 1
delta        : [0.8 0.8 0.8]
Mdc          : [[-1.  0.  0.]
 [ 0.  0. -1.]
 [ 0.  1.  0.]]
Pxyz_c       : [0.40000153 0.40000153 0.40000153]
tr           : 1000.0
flip_angle   : 0.0
te           : 0.0
ti           : 0.0
fov          : 256.0

In a GIFTI surface file, it appears as follows:

>>> dict(gii.darrays[0].meta)
{[...]
 'VolGeomWidth': '320',
 'VolGeomHeight': '320',
 'VolGeomDepth': '320',
 'VolGeomXsize': '0.800000',
 'VolGeomYsize': '0.800000',
 'VolGeomZsize': '0.800000',
 'VolGeomX_R': '-1.000000',
 'VolGeomX_A': '0.000000',
 'VolGeomX_S': '0.000000',
 'VolGeomY_R': '0.000000',
 'VolGeomY_A': '0.000000',
 'VolGeomY_S': '-1.000000',
 'VolGeomZ_R': '0.000000',
 'VolGeomZ_A': '1.000000',
 'VolGeomZ_S': '0.000000',
 'VolGeomC_R': '0.400002',
 'VolGeomC_A': '0.400002',
 'VolGeomC_S': '0.400002',
 [...]
}

Equivalent structures also appear in the footer of surface files:

>>> lh_white = nb.freesurfer.read_geometry('surf/lh.white', read_metadata=True)
>>> lh_white[2]
OrderedDict([('head', array([ 2,  0, 20], dtype=int32)),
             ('valid', '1  # volume info valid'),
             ('filename', '../mri/filled-pretess255.mgz'),
             ('volume', array([320, 320, 320])),
             ('voxelsize', array([0.80000001, 0.80000001, 0.80000001])),
             ('xras', array([-1.,  0.,  0.])),
             ('yras', array([ 0.,  0., -1.])),
             ('zras', array([0., 1., 0.])),
             ('cras', array([0.40000153, 0.40000153, 0.40000153]))])

The LTA affine transform file format also contains a representation of the volume geometry, but I won't go into that here.

Here are functions to create VolGeoms from each structure:

from operator import itemgetter
import numpy as np

class VolGeom:
    ...

    @classmethod
    def from_mghheader(cls, header):
        return cls(
            shape=header.dims[:3],
            zooms=header.delta,
            cosines=header.Mdc.T,
            c_ras=header.Pxyz_c,
        )

    @classmethod
    def from_surf_footer(cls, footer):
        return cls(
            shape=footer['volume],
            zooms=footer['voxelsize'],
            cosines=np.vstack(itemgetter('xras', 'yras', 'zras')(footer)).T,
            c_ras=footer['cras'],
        )

    @classmethod
    def from_gifti_metadata(cls, meta):
        return cls(
            shape=np.int16((meta['VolGeomWidth'], meta['VolGeomHeight'], meta['VolGeomDepth'])),
            zooms=np.float64((meta['VolGeomXsize'], meta['VolGeomYsize'], meta['VolGeomZsize'])),
            cosines=np.float64(
                itemgetter(*(f'VolGeom{col}_{row}' for row in 'RAS' for col in 'XYZ'))(meta)
            ).reshape((3, 3)),
            c_ras=np.float64(meta['VolGeomC_R'], meta['VolGeomC_A'], meta['VolGeomC_S']),
        )

Affines

We can define the tkregister and scanner affines as follows:

import nibabel as nb

class VolGeom:
    ...

    def tkreg_affine(self):
        tkrcosines = np.array([[-1, 0, 0], [0, 0, 1], [0, -1, 0]])
        mat = tkrcosines * self.zooms
        return nb.affines.from_matvec(mat, -mat @ self.shape / 2)

    def scanner_affine(self):
        mat = self.cosines * self.zooms
        return nb.affines.from_matvec(mat, self.c_ras - mat @ self.shape / 2)

Note that these would be the same if c_ras == [0, 0, 0] and self.cosines == tkrcosines. When the image is non-oblique, the difference is the offset of the image center from the origin (in mm), or c_ras. Some tools (HCP pipelines, older versions of fMRIPrep) have mistakenly added c_ras to TKR coordinates, which works for non-oblique datasets.

To find the transform from RAS coordinates in tkregister space to RAS coordinates in scanner space, we can multiply the scanner affine by the inverse of the tkregister affine:

class VolGeom:
    ...

    def tkreg2scanner(self):
        return self.scanner_affine() @ np.linalg.inv(self.tkreg_affine())

Full definition

from dataclasses import dataclass
from operator import itemgetter

import numpy as np

import nibabel as nb

@dataclass
class VolGeom:
    # Using tuples to show shape. Will mostly be numpy arrays.
    shape: tuple[int, int, int]
    zooms: tuple[float, float, float]
    cosines: tuple[tuple[float, float, float], tuple[float, float, float], tuple[float, float, float]]
    c_ras: tuple[float, float, float]

    @classmethod
    def from_mghheader(cls, header):
        return cls(
            shape=header.dims[:3],
            zooms=header.delta,
            cosines=header.Mdc.T,
            c_ras=header.Pxyz_c,
        )

    @classmethod
    def from_surf_footer(cls, footer):
        return cls(
            shape=footer['volume'],
            zooms=footer['voxelsize'],
            cosines=np.vstack(itemgetter('xras', 'yras', 'zras')(footer)).T,
            c_ras=footer['cras'],
        )

    @classmethod
    def from_gifti_metadata(cls, meta):
        return cls(
            shape=np.int16((meta['VolGeomWidth'], meta['VolGeomHeight'], meta['VolGeomDepth'])),
            zooms=np.float64((meta['VolGeomXsize'], meta['VolGeomYsize'], meta['VolGeomZsize'])),
            cosines=np.float64(
                itemgetter(*(f'VolGeom{col}_{row}' for row in 'RAS' for col in 'XYZ'))(meta)
            ).reshape((3, 3)),
            c_ras=np.float64((meta['VolGeomC_R'], meta['VolGeomC_A'], meta['VolGeomC_S'])),
    )

    def tkreg_affine(self):
        tkrcosines = np.array([[-1, 0, 0], [0, 0, 1], [0, -1, 0]])
        mat = tkrcosines * self.zooms
        return nb.affines.from_matvec(mat, -mat @ self.shape / 2)

    def scanner_affine(self):
        mat = self.cosines * self.zooms
        return nb.affines.from_matvec(mat, self.c_ras - mat @ self.shape / 2)

    def tkreg2scanner(self):
        return self.scanner_affine() @ np.linalg.inv(self.tkreg_affine())

Converting a FreeSurfer GIFTI surface to RAS coordinates

Convert lh.white to GIFTI with and without --to-scanner:

$ mris_convert surf/lh.white hemi-L_white.surf.gii
$ mris_convert --to-scanner surf/lh.white hemi-L_desc-scanner_white.surf.gii

Now we can perform our own --to-scanner conversion:

orig = nb.load('hemi-L_white.surf.gii')
coords, triangles = orig.darrays
geom = VolGeom.from_gifti_metadata(coords.meta)

scanner_coords = coords.__class__(
    nb.affines.apply_affine(geom.tkreg2scanner(), coords.data).astype('f4'),
    meta=coords.meta
)
ourscanner = nb.GiftiImage(darrays=[scanner_coords, triangles], meta=orig.meta)

And compare to FreeSurfer's:

scanner = nb.load('hemi-L_desc-scanner_white.surf.gii')
assert np.allclose(our_scanner.darrays[0].data, scanner.darrays[0].data, atol=1e-4)