07. Save and load T1-weighted MRI scan along with anatomical landmarks in BIDS

When working with MEEG data in the domain of source localization, we usually have to deal with aligning several coordinate systems, such as the coordinate systems of …

  • the head of a study participant

  • the recording device (in the case of MEG)

  • the anatomical MRI scan of a study participant

The process of aligning these frames is also called coregistration, and is performed with the help of a transformation matrix, called trans in MNE.

In this tutorial, we show how MNE-BIDS can be used to save a T1 weighted MRI scan in BIDS format, and to encode all information of the trans object in a BIDS compatible way.

Finally, we will automatically reproduce our trans object from a BIDS directory.

See the documentation pages in the MNE docs for more information on source alignment and coordinate frames

Note

For this example you will need to install matplotlib and nilearn on top of your usual mne-bids installation.

# Authors: Stefan Appelhoff <stefan.appelhoff@mailbox.org>
#          Alex Rockhill <aprockhill206@gmail.com>
#          Alex Gramfort <alexandre.gramfort@inria.fr>
#
# License: BSD (3-clause)

Let’s import everything we need for this example:

import os.path as op
import shutil

import numpy as np
import matplotlib.pyplot as plt

from nilearn.plotting import plot_anat

import mne
from mne.datasets import sample
from mne import head_to_mri

from mne_bids import (write_raw_bids, BIDSPath, write_anat, get_anat_landmarks,
                      get_head_mri_trans, print_dir_tree)

Out:

/home/stefanappelhoff/miniconda3/envs/mne_bids/lib/python3.8/site-packages/nilearn/datasets/__init__.py:87: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
  warn("Fetchers from the nilearn.datasets module will be "

We will be using the MNE sample data and write a basic BIDS dataset. For more information, you can checkout the respective example.

data_path = sample.data_path()
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2, 'Visual/Left': 3,
            'Visual/Right': 4, 'Smiley': 5, 'Button': 32}
raw_fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
events_data = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw-eve.fif')
output_path = op.abspath(op.join(data_path, '..', 'MNE-sample-data-bids'))
fs_subjects_dir = op.join(data_path, 'subjects')  # FreeSurfer subjects dir

To ensure the output path doesn’t contain any leftover files from previous tests and example runs, we simply delete it.

Warning

Do not delete directories that may contain important data!

Read the input data and store it as BIDS data.

raw = mne.io.read_raw_fif(raw_fname)
raw.info['line_freq'] = 60  # specify power line frequency as required by BIDS

sub = '01'
ses = '01'
task = 'audiovisual'
run = '01'
bids_path = BIDSPath(subject=sub, session=ses, task=task,
                     run=run, root=output_path)
write_raw_bids(raw, bids_path, events_data=events_data,
               event_id=event_id, overwrite=True)

Out:

Opening raw data file /home/stefanappelhoff/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.
Opening raw data file /home/stefanappelhoff/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/README'...

References
----------
Appelhoff, S., Sanderson, M., Brooks, T., Vliet, M., Quentin, R., Holdgraf, C., Chaumon, M., Mikulan, E., Tavabi, K., Höchenberger, R., Welke, D., Brunner, C., Rockhill, A., Larson, E., Gramfort, A. and Jas, M. (2019). MNE-BIDS: Organizing electrophysiological data into the BIDS format and facilitating their analysis. Journal of Open Source Software 4: (1896). https://doi.org/10.21105/joss.01896

Niso, G., Gorgolewski, K. J., Bock, E., Brooks, T. L., Flandin, G., Gramfort, A., Henson, R. N., Jas, M., Litvak, V., Moreau, J., Oostenveld, R., Schoffelen, J., Tadel, F., Wexler, J., Baillet, S. (2018). MEG-BIDS, the brain imaging data structure extended to magnetoencephalography. Scientific Data, 5, 180110. https://doi.org/10.1038/sdata.2018.110


Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/participants.tsv'...

participant_id  age     sex     hand
sub-01  n/a     n/a     n/a

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/participants.json'...

{
    "participant_id": {
        "Description": "Unique participant identifier"
    },
    "age": {
        "Description": "Age of the participant at time of testing",
        "Units": "years"
    },
    "sex": {
        "Description": "Biological sex of the participant",
        "Levels": {
            "F": "female",
            "M": "male"
        }
    },
    "hand": {
        "Description": "Handedness of the participant",
        "Levels": {
            "R": "right",
            "L": "left",
            "A": "ambidextrous"
        }
    }
}

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_coordsystem.json'...

{
    "MEGCoordinateSystem": "ElektaNeuromag",
    "MEGCoordinateUnits": "m",
    "MEGCoordinateSystemDescription": "RAS orientation and the origin between the ears",
    "HeadCoilCoordinates": {
        "NAS": [
            3.725290298461914e-09,
            0.10260561108589172,
            4.190951585769653e-09
        ],
        "LPA": [
            -0.07137660682201385,
            0.0,
            5.122274160385132e-09
        ],
        "RPA": [
            0.07526767998933792,
            0.0,
            5.587935447692871e-09
        ],
        "coil1": [
            0.032922741025686264,
            0.09897983074188232,
            0.07984329760074615
        ],
        "coil2": [
            -0.06998106092214584,
            0.06771647930145264,
            0.06888450682163239
        ],
        "coil3": [
            -0.07260829955339432,
            -0.02086828649044037,
            0.0971473976969719
        ],
        "coil4": [
            0.04996863007545471,
            -0.007233052980154753,
            0.1228904277086258
        ]
    },
    "HeadCoilCoordinateSystem": "ElektaNeuromag",
    "HeadCoilCoordinateUnits": "m",
    "AnatomicalLandmarkCoordinates": {
        "NAS": [
            3.725290298461914e-09,
            0.10260561108589172,
            4.190951585769653e-09
        ],
        "LPA": [
            -0.07137660682201385,
            0.0,
            5.122274160385132e-09
        ],
        "RPA": [
            0.07526767998933792,
            0.0,
            5.587935447692871e-09
        ]
    },
    "AnatomicalLandmarkCoordinateSystem": "ElektaNeuromag",
    "AnatomicalLandmarkCoordinateUnits": "m"
}

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_events.tsv'...

onset   duration        trial_type      value   sample
3.6246181587150867      0.0     Auditory/Right  2       2177
4.237323479067476       0.0     Visual/Left     3       2545
4.946596485779753       0.0     Auditory/Left   1       2971
5.692498614904401       0.0     Visual/Right    4       3419
6.41342634238425        0.0     Auditory/Right  2       3852

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/dataset_description.json'...

{
    "Name": " ",
    "BIDSVersion": "1.6.0",
    "DatasetType": "raw",
    "Authors": [
        "[Unspecified]"
    ]
}

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_meg.json'...

{
    "TaskName": "audiovisual",
    "Manufacturer": "Elekta",
    "PowerLineFrequency": 60,
    "SamplingFrequency": 600.614990234375,
    "SoftwareFilters": "n/a",
    "RecordingDuration": 277.7136813300495,
    "RecordingType": "continuous",
    "DewarPosition": "n/a",
    "DigitizedLandmarks": true,
    "DigitizedHeadPoints": true,
    "MEGChannelCount": 306,
    "MEGREFChannelCount": 0,
    "ContinuousHeadLocalization": false,
    "HeadCoilFrequency": [],
    "EEGChannelCount": 60,
    "EOGChannelCount": 1,
    "ECGChannelCount": 0,
    "EMGChannelCount": 0,
    "MiscChannelCount": 0,
    "TriggerChannelCount": 9
}

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_channels.tsv'...

name    type    units   low_cutoff      high_cutoff     description     sampling_frequency      status  status_description
MEG 0113        MEGGRADPLANAR   T/m     0.10000000149011612     172.17630004882812      Planar Gradiometer      600.614990234375        good    n/a
MEG 0112        MEGGRADPLANAR   T/m     0.10000000149011612     172.17630004882812      Planar Gradiometer      600.614990234375        good    n/a
MEG 0111        MEGMAG  T       0.10000000149011612     172.17630004882812      Magnetometer    600.614990234375        good    n/a
MEG 0122        MEGGRADPLANAR   T/m     0.10000000149011612     172.17630004882812      Planar Gradiometer      600.614990234375        good    n/a
MEG 0123        MEGGRADPLANAR   T/m     0.10000000149011612     172.17630004882812      Planar Gradiometer      600.614990234375        good    n/a
Copying data files to sub-01_ses-01_task-audiovisual_run-01_meg.fif
Reserving possible split file sub-01_ses-01_task-audiovisual_run-01_split-01_meg.fif
Writing /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_meg.fif
Closing /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_meg.fif
[done]

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/sub-01_ses-01_scans.tsv'...

filename        acq_time
meg/sub-01_ses-01_task-audiovisual_run-01_meg.fif       2002-12-03T19:01:10.720100Z
Wrote /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/sub-01_ses-01_scans.tsv entry with meg/sub-01_ses-01_task-audiovisual_run-01_meg.fif.

BIDSPath(
root: /home/stefanappelhoff/mne_data/MNE-sample-data-bids
datatype: meg
basename: sub-01_ses-01_task-audiovisual_run-01_meg.fif)

Print the directory tree

Out:

|MNE-sample-data-bids/
|--- README
|--- dataset_description.json
|--- participants.json
|--- participants.tsv
|--- sub-01/
|------ ses-01/
|--------- sub-01_ses-01_scans.tsv
|--------- meg/
|------------ sub-01_ses-01_coordsystem.json
|------------ sub-01_ses-01_task-audiovisual_run-01_channels.tsv
|------------ sub-01_ses-01_task-audiovisual_run-01_events.tsv
|------------ sub-01_ses-01_task-audiovisual_run-01_meg.fif
|------------ sub-01_ses-01_task-audiovisual_run-01_meg.json

Writing T1 image

Now let’s assume that we have also collected some T1 weighted MRI data for our subject. And furthermore, that we have already aligned our coordinate frames (using e.g., the coregistration GUI) and obtained a transformation matrix trans.

# Get the path to our MRI scan
t1_fname = op.join(fs_subjects_dir, 'sample', 'mri', 'T1.mgz')

# Load the transformation matrix and show what it looks like
trans_fname = op.join(data_path, 'MEG', 'sample',
                      'sample_audvis_raw-trans.fif')
trans = mne.read_trans(trans_fname)
print(trans)

Out:

<Transform | head->MRI (surface RAS)>
[[ 0.99930954  0.01275934  0.0348942   0.00206991]
 [ 0.00998479  0.81240475 -0.58300853  0.01130214]
 [-0.03578702  0.58295429  0.81171638 -0.02755522]
 [ 0.          0.          0.          1.        ]]

We can save the MRI to our existing BIDS directory and at the same time create a JSON sidecar file that contains metadata, we will later use to retrieve our transformation matrix trans. The metadata will here consist of the coordinates of three anatomical landmarks (LPA, Nasion and RPA (=left and right preauricular points) expressed in voxel coordinates w.r.t. the T1 image.

# First create the BIDSPath object.
t1w_bids_path = BIDSPath(subject=sub, session=ses, root=output_path,
                         suffix='T1w')

# use ``trans`` to transform landmarks from the ``raw`` file to
# the voxel space of the image
landmarks = get_anat_landmarks(
    t1_fname,  # path to the MRI scan
    info=raw.info,  # the MEG data file info from the same subject as the MRI
    trans=trans,  # our transformation matrix
    fs_subject='sample',  # FreeSurfer subject
    fs_subjects_dir=fs_subjects_dir,  # FreeSurfer subjects directory
)

# We use the write_anat function
t1w_bids_path = write_anat(
    image=t1_fname,  # path to the MRI scan
    bids_path=t1w_bids_path,
    landmarks=landmarks,  # the landmarks in MRI voxel space
    verbose=True  # this will print out the sidecar file
)
anat_dir = t1w_bids_path.directory

Out:

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/anat/sub-01_ses-01_T1w.json'...

{
    "AnatomicalLandmarkCoordinates": {
        "LPA": [
            197.25741767555172,
            153.00085676768788,
            138.58946456240983
        ],
        "NAS": [
            124.62090970591665,
            95.74083306302882,
            222.65943140321255
        ],
        "RPA": [
            50.71438289229903,
            158.2488189432002,
            140.0536763406783
        ]
    }
}

Let’s have another look at our BIDS directory

Out:

|MNE-sample-data-bids/
|--- README
|--- dataset_description.json
|--- participants.json
|--- participants.tsv
|--- sub-01/
|------ ses-01/
|--------- sub-01_ses-01_scans.tsv
|--------- anat/
|------------ sub-01_ses-01_T1w.json
|------------ sub-01_ses-01_T1w.nii.gz
|--------- meg/
|------------ sub-01_ses-01_coordsystem.json
|------------ sub-01_ses-01_task-audiovisual_run-01_channels.tsv
|------------ sub-01_ses-01_task-audiovisual_run-01_events.tsv
|------------ sub-01_ses-01_task-audiovisual_run-01_meg.fif
|------------ sub-01_ses-01_task-audiovisual_run-01_meg.json

Our BIDS dataset is now ready to be shared. We can easily estimate the transformation matrix using MNE-BIDS and the BIDS dataset. This function converts the anatomical landmarks stored in the T1 sidecar file into FreeSurfer surface RAS space, and aligns the landmarks in the electrophysiology data with them. This way your electrophysiology channel locations can be transformed to surface RAS space using the trans which is crucial for source localization and other uses of the FreeSurfer surfaces.

Note

If this dataset were shared with you, you would first have to use the T1 image as input for the FreeSurfer recon-all, see How MNE uses FreeSurfer’s outputs.

Out:

Opening raw data file /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_meg.fif...
    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
    Range : 25800 ... 192599 =     42.956 ...   320.670 secs
Ready.
Reading events from /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_events.tsv.
Reading channel info from /home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/meg/sub-01_ses-01_task-audiovisual_run-01_channels.tsv.
/home/stefanappelhoff/Desktop/bids/mne-bids/mne_bids/read.py:529: RuntimeWarning: The unit for channel(s) STI 001, STI 002, STI 003, STI 004, STI 005, STI 006, STI 014, STI 015, STI 016 has changed from V to NA.
  raw.set_channel_types(ch_type_map_avail)

Finally, let’s use the T1 weighted MRI image and plot the anatomical landmarks Nasion, LPA, and RPA onto the brain image. For that, we can extract the location of Nasion, LPA, and RPA from the MEG file, apply our transformation matrix trans, and plot the results.

# Get Landmarks from MEG file, 0, 1, and 2 correspond to LPA, NAS, RPA
# and the 'r' key will provide us with the xyz coordinates. The coordinates
# are expressed here in MEG Head coordinate system.
pos = np.asarray((raw.info['dig'][0]['r'],
                  raw.info['dig'][1]['r'],
                  raw.info['dig'][2]['r']))

# We now use the ``head_to_mri`` function from MNE-Python to convert MEG
# coordinates to MRI scanner RAS space. For the conversion we use our
# estimated transformation matrix and the MEG coordinates extracted from the
# raw file. `subjects` and `subjects_dir` are used internally, to point to
# the T1-weighted MRI file: `t1_mgh_fname`. Coordinates are is mm.
mri_pos = head_to_mri(pos=pos,
                      subject='sample',
                      mri_head_t=estim_trans,
                      subjects_dir=fs_subjects_dir)

# Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function
t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')

# Plot it
fig, axs = plt.subplots(3, 1, figsize=(7, 7), facecolor='k')
for point_idx, label in enumerate(('LPA', 'NAS', 'RPA')):
    plot_anat(t1_nii_fname, axes=axs[point_idx],
              cut_coords=mri_pos[point_idx, :],
              title=label, vmax=160)
plt.show()
convert mri and trans

Writing FLASH MRI image

We can write another types of MRI data such as FLASH images for BEM models

flash_fname = op.join(fs_subjects_dir, 'sample', 'mri', 'flash', 'mef05.mgz')

flash_bids_path = \
    BIDSPath(subject=sub, session=ses, root=output_path, suffix='FLASH')

write_anat(
    image=flash_fname,
    bids_path=flash_bids_path,
    verbose=True
)

Out:

BIDSPath(
root: /home/stefanappelhoff/mne_data/MNE-sample-data-bids
datatype: anat
basename: sub-01_ses-01_FLASH.nii.gz)

Writing defaced and anonymized T1 image

We can deface the MRI for anonymization by passing deface=True.

t1w_bids_path = write_anat(
    image=t1_fname,  # path to the MRI scan
    bids_path=bids_path,
    landmarks=landmarks,
    deface=True,
    overwrite=True,
    verbose=True  # this will print out the sidecar file
)
anat_dir = t1w_bids_path.directory

# Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function
t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')

# Plot it
fig, ax = plt.subplots()
plot_anat(t1_nii_fname, axes=ax, title='Defaced', vmax=160)
plt.show()
convert mri and trans

Out:

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/anat/sub-01_ses-01_T1w.json'...

{
    "AnatomicalLandmarkCoordinates": {
        "LPA": [
            197.25741767555172,
            153.00085676768788,
            138.58946456240983
        ],
        "NAS": [
            124.62090970591665,
            95.74083306302882,
            222.65943140321255
        ],
        "RPA": [
            50.71438289229903,
            158.2488189432002,
            140.0536763406783
        ]
    }
}

Writing defaced and anonymized FLASH MRI image

Defacing the FLASH works just like the T1 as long as they are aligned.

# use ``trans`` to transform landmarks from the ``raw`` file to
# the voxel space of the image
landmarks = get_anat_landmarks(
    flash_fname,  # path to the FLASH scan
    info=raw.info,  # the MEG data file info from the same subject as the MRI
    trans=trans,  # our transformation matrix
    fs_subject='sample',  # freesurfer subject
    fs_subjects_dir=fs_subjects_dir,  # freesurfer subjects directory
)

flash_bids_path = write_anat(
    image=flash_fname,  # path to the MRI scan
    bids_path=flash_bids_path,
    landmarks=landmarks,
    deface=True,
    overwrite=True,
    verbose=True  # this will print out the sidecar file
)

# Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function
flash_nii_fname = op.join(anat_dir, 'sub-01_ses-01_FLASH.nii.gz')

# Plot it
fig, ax = plt.subplots()
plot_anat(flash_nii_fname, axes=ax, title='Defaced', vmax=700)
plt.show()
convert mri and trans

Out:

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/anat/sub-01_ses-01_FLASH.json'...

{
    "AnatomicalLandmarkCoordinates": {
        "LPA": [
            124.51988917211392,
            156.7554808781791,
            164.5250043352138
        ],
        "NAS": [
            28.620126675736444,
            116.1185261240982,
            95.42797291038784
        ],
        "RPA": [
            120.56940719564132,
            157.75593669628364,
            17.937353558892525
        ]
    }
}

Using manual landmark coordinates in scanner RAS

You can also find landmarks with a 3D image viewer (e.g. FreeView) if you have not aligned the channel locations (including fiducials) using the coregistration GUI or if this is just more convenient.

Note

In FreeView, you need to use “RAS” and not “TkReg RAS” for this. You can also use voxel coordinates but, in FreeView, they are integers and so not as precise as the “RAS” decimal numbers.

flash_ras_landmarks = \
    np.array([[-74.53102838, 19.62854953, -52.2888194],
              [-1.89454315, 103.69850925, 4.97120376],
              [72.01200673, 21.09274883, -57.53678375]]) / 1e3  # mm -> m

landmarks = mne.channels.make_dig_montage(
    lpa=flash_ras_landmarks[0],
    nasion=flash_ras_landmarks[1],
    rpa=flash_ras_landmarks[2],
    coord_frame='ras'
)

flash_bids_path = write_anat(
    image=flash_fname,  # path to the MRI scan
    bids_path=flash_bids_path,
    landmarks=landmarks,
    deface=True,
    overwrite=True,
    verbose=True  # this will print out the sidecar file
)

# Plot it
fig, ax = plt.subplots()
plot_anat(flash_nii_fname, axes=ax, title='Defaced', vmax=700)
plt.show()
convert mri and trans

Out:

Writing '/home/stefanappelhoff/mne_data/MNE-sample-data-bids/sub-01/ses-01/anat/sub-01_ses-01_FLASH.json'...

{
    "AnatomicalLandmarkCoordinates": {
        "LPA": [
            124.51988917137612,
            156.75548087980883,
            164.525004336295
        ],
        "NAS": [
            28.62012667241129,
            116.11852612679817,
            95.42797290948731
        ],
        "RPA": [
            120.56940719949279,
            157.75593669708275,
            17.937353558283036
        ]
    }
}

Total running time of the script: ( 0 minutes 33.259 seconds)

Gallery generated by Sphinx-Gallery