How MNE uses FreeSurfer's outputs

This tutorial explains how MRI coordinate frames are handled in MNE-Python, and how MNE-Python integrates with FreeSurfer for handling MRI data and source space data in general.

As usual we'll start by importing the necessary packages; for this tutorial that includes :mod:nibabel to handle loading the MRI images (MNE-Python also uses :mod:nibabel under the hood). We'll also use a special :mod:`Matplotlib

<matplotlib.patheffects>` function for adding outlines to text, so that text is readable on top of an MRI image.

MRI coordinate frames

Let's start out by looking at the sample subject MRI. Following standard FreeSurfer convention, we look at :file:T1.mgz, which gets created from the original MRI :file:sample/mri/orig/001.mgz when you run the FreeSurfer command recon-all <https://surfer.nmr.mgh.harvard.edu/fswiki/recon-all>_. Here we use :mod:nibabel to load the T1 image, and the resulting object's :meth:~nibabel.spatialimages.SpatialImage.orthoview method to view it.

Notice that the axes in the :meth:~nibabel.spatialimages.SpatialImage.orthoview figure are labeled L-R, S-I, and P-A. These reflect the standard RAS (right-anterior-superior) coordinate system that is widely used in MRI imaging. If you are unfamiliar with RAS coordinates, see the excellent nibabel tutorial :doc:nibabel:coordinate_systems.

Nibabel already takes care of some coordinate frame transformations under the hood, so let's do it manually so we understand what is happening. First let's get our data as a 3D array and note that it's already a standard size:

These data are voxel intensity values. Here they are unsigned integers in the range 0-255, though in general they can be floating point values. A value data[i, j, k] at a given index triplet (i, j, k) corresponds to some real-world physical location (x, y, z) in space. To get its physical location, first we have to choose what coordinate frame we're going to use.

For example, we could choose a geographical coordinate frame, with origin is at the center of the earth, Z axis through the north pole, X axis through the prime meridian (zero degrees longitude), and Y axis orthogonal to these forming a right-handed coordinate system. This would not be a very useful choice for defining the physical locations of the voxels during the MRI acquisition for analysis, but you could nonetheless figure out the transformation that related the (i, j, k) to this coordinate frame.

Instead, each scanner defines a more practical, native coordinate system that it uses during acquisition, usually related to the physical orientation of the scanner itself and/or the subject within it. During acquisition the relationship between the voxel indices (i, j, k) and the physical location (x, y, z) in the scanner's native coordinate frame is saved in the image's affine transformation.

.. sidebar:: Under the hood

``mne.transforms.apply_trans`` effectively does a matrix multiplication
(i.e., :func:`numpy.dot`), with a little extra work to handle the shape
mismatch (the affine has shape ``(4, 4)`` because it includes a
*translation*, which is applied separately).

We can use :mod:nibabel to examine this transformation, keeping in mind that it processes everything in units of millimeters, unlike MNE where things are always in SI units (meters).

This allows us to take an arbitrary voxel or slice of data and know where it is in the scanner's native physical space (x, y, z) (in mm) by applying the affine transformation to the voxel coordinates.

If you have a point (x, y, z) in scanner-native RAS space and you want the corresponding voxel number, you can get it using the inverse of the affine. This involves some rounding, so it's possible to end up off by one voxel if you're not careful:

Let's write a short function to visualize where our voxel lies in an image, and annotate it in RAS space (rounded to the nearest millimeter):

Notice that the axis scales (i, j, and k) are still in voxels (ranging from 0-255); it's only the annotation text that we've translated into real-world RAS in millimeters.

"MRI coordinates" in MNE-Python: FreeSurfer surface RAS

While :mod:nibabel uses scanner RAS (x, y, z) coordinates, FreeSurfer uses a slightly different coordinate frame: MRI surface RAS. The transform from voxels to the FreeSurfer MRI surface RAS coordinate frame is known in the FreeSurfer documentation <https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems>_ as Torig, and in nibabel as :meth:`vox2ras_tkr

<nibabel.freesurfer.mghformat.MGHHeader.get_vox2ras_tkr>. This transformation sets the center of its coordinate frame in the middle of the conformed volume dimensions (``N / 2.``) with the axes oriented along the axes of the volume itself. For more information, seecoordinate_systems`.

Since MNE-Python uses FreeSurfer extensively for surface computations (e.g., white matter, inner/outer skull meshes), internally MNE-Python uses the Freeurfer surface RAS coordinate system (not the :mod:nibabel scanner RAS system) for as many computations as possible, such as all source space and BEM mesh vertex definitions.

Whenever you see "MRI coordinates" or "MRI coords" in MNE-Python's documentation, you should assume that we are talking about the "FreeSurfer MRI surface RAS" coordinate frame!

We can do similar computations as before to convert the given voxel indices into FreeSurfer MRI coordinates (i.e., what we call "MRI coordinates" or "surface RAS" everywhere else in MNE), just like we did above to convert voxel indices to scanner RAS:

Knowing these relationships and being mindful about transformations, we can get from a point in any given space to any other space. Let's start out by plotting the Nasion on a saggital MRI slice:

When we print the nasion, it displays as a DigPoint and shows its coordinates in millimeters, but beware that the underlying data is actually stored in meters <units>, so before transforming and plotting we'll convert to millimeters:

We can also take the digitization point from the MEG data, which is in the "head" coordinate frame.

Let's look at the nasion in the head coordinate frame:

.. sidebar:: Head coordinate frame

 The head coordinate frame in MNE is the "Neuromag" head coordinate
 frame. The origin is given by the intersection between a line connecting
 the LPA and RPA and the line orthogonal to it that runs through the
 nasion. It is also in RAS orientation, meaning that +X runs through
 the RPA, +Y goes through the nasion, and +Z is orthogonal to these
 pointing upward. See `coordinate_systems` for more information.

Notice that in "head" coordinate frame the nasion has values of 0 for the x and z directions (which makes sense given that the nasion is used to define the y axis in that system). To convert from head coordinate frame to voxels, we first apply the head → MRI (surface RAS) transform from a :file:trans file (typically created with the MNE-Python coregistration GUI), then convert meters → millimeters, and finally apply the inverse of Torig to get to voxels.

Under the hood, functions like :func:mne.setup_source_space, :func:mne.setup_volume_source_space, and :func:mne.compute_source_morph make extensive use of these coordinate frames.

Using FreeSurfer's surface reconstructions

An important part of what FreeSurfer does is provide cortical surface reconstructions. For example, let's load and view the white surface of the brain. This is a 3D mesh defined by a set of vertices (conventionally called rr) with shape (n_vertices, 3) and a set of triangles (tris) with shape (n_tris, 3) defining which vertices in rr form each triangular facet of the mesh.

Let's actually plot it:

We can also plot the mesh on top of an MRI slice. The mesh surfaces are defined in millimeters in the MRI (FreeSurfer surface RAS) coordinate frame, so we can convert them to voxels by applying the inverse of the Torig transform:

This is the method used by :func:mne.viz.plot_bem to show the BEM surfaces.

Cortical alignment (spherical)

A critical function provided by FreeSurfer is spherical surface alignment of cortical surfaces, maximizing sulcal-gyral alignment. FreeSurfer first expands the cortical surface to a sphere, then aligns it optimally with fsaverage. Because the vertex ordering is preserved when expanding to a sphere, a given vertex in the source (sample) mesh can be mapped easily to the same location in the destination (fsaverage) mesh, and vice-versa.

Let's look a bit more closely at the spherical alignment by overlaying the two spherical meshes as wireframes and zooming way in (the purple points are separated by about 1 mm):

You can see that the fsaverage (purple) mesh is uniformly spaced, and the mesh for subject "sample" (in cyan) has been deformed along the spherical surface by FreeSurfer. This deformation is designed to optimize the sulcal-gyral alignment.

Surface decimation

These surfaces have a lot of vertices, and in general we only need to use a subset of these vertices for creating source spaces. A uniform sampling can easily be achieved by subsampling in the spherical space. To do this, we use a recursively subdivided icosahedron or octahedron. For example, let's load a standard oct-6 source space, and at the same zoom level as before visualize how it subsampled the dense mesh:

We can also then look at how these two meshes compare by plotting the original, high-density mesh as well as our decimated mesh white surfaces.