Obtain head -> MRI transform for the MOUS dataset (DSC_3011020.09_236) for use with MNE-Python.

mous_trans.py

"""
Obtain head -> MRI transform for the MOUS dataset for use with MNE-Python.

If you want to analyze the "Mother of all unification studies" dataset with
FreeSurfer and MNE-Python, you will need the MRI->HEAD coordinate
transformations. In the original dataset, the T1w_space_CTF.nii files contain
this information, but it takes some doing to distill this into a transform
object as used by MNE-Python.

This script assumes you have ran FreeSurfer on the T1w.nii files of all the
subjects.

Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
"""
import argparse

import mne
import nibabel
import numpy as np
from mne.transforms import Transform

# Set these paths to match your setup
#####################################

# Path to where you've downloaded the MOUS data to
mous_path = "/path/to/DSC_3011020.09_236"

# Path to your FreeSurfer SUBJECTS_DIR
subjects_dir = "/path/to/DSC_3011020.09_236/derivatives/freesurfer"

# Subject to compute the transform for
subject = "sub-V1001"

# Original T1w.nii file
fname_anat = f"{mous_path}/{subject}/anat/T1w.nii"

# Provided T1w_space_CTF.nii file
fname_anat_ctf = f"{mous_path}/{subject}/anat/T1w_space_CTF.nii"

# FreeSurfer T1 file
fname_freesurfer_t1 = f"{subjects_dir}/{subject}/mri/T1.mgz"

# One of the files containing raw MEG data
fname_ctf = f"{mous_path}/{subject}/meg/{subject}_task-rest_meg.ds"

# Where to save the trans.fif file to
fname_trans = f"{subjects_dir}/{subject}/bem/{subject}-trans.fif"


# Do coordinate transforms until we get what we want.
#####################################################

# Freesurfer RAS -> Freesurfer voxel
freesurfer_mri = nibabel.load(fname_freesurfer_t1)
freesurfer_ras_to_freesurfer_voxel = Transform(
    "mri", "mri_voxel", np.linalg.inv(freesurfer_mri.header.get_vox2ras_tkr())
)

# Freesurfer voxel (T1.mgz) -> BIDS voxel (sub-????_T1w.nii) transform
bids_mri = nibabel.load(fname_anat)
freesurfer_voxel_to_scanner_ras = Transform("mri_voxel", "mri", freesurfer_mri.affine)
scanner_ras_to_bids_voxel = Transform(
    "mri", "mri_voxel", np.linalg.inv(bids_mri.affine)
)

# BIDS voxel -> CTF voxel transform
bids_voxel_to_ctf_voxel = Transform(
    "mri_voxel",
    "mri_voxel",
    np.array(
        [
            [0, -1, 0, 255],
            [0, 0, 1, 0],
            [1, 0, 0, 0],
            [0, 0, 0, 1],
        ]
    ),
)

# CTF voxel to CTF space transform
ctf_mri = nibabel.load(fname_anat_ctf)
ctf_voxel_to_ctf = mne.transforms.Transform("mri_voxel", "ctf_meg", ctf_mri.affine)

# CTF space to head space transform
info = mne.io.read_raw_ctf(fname_ctf).info
ctf_to_head = info["ctf_head_t"]

# Put together the final transform
transforms = [
    freesurfer_ras_to_freesurfer_voxel,
    freesurfer_voxel_to_scanner_ras,
    scanner_ras_to_bids_voxel,
    bids_voxel_to_ctf_voxel,
    ctf_voxel_to_ctf,
    ctf_to_head,
]

mri_to_head = np.eye(4)
for t in transforms[::-1]:
    mri_to_head = mri_to_head @ t["trans"]
mri_to_head[:3, 3] /= 1000
head_to_mri = Transform("head", "mri", np.linalg.inv(mri_to_head))
mne.write_trans(fname_trans, head_to_mri, overwrite=True)

fig = mne.viz.plot_alignment(
    info=info,
    trans=head_to_mri,
    subject=subject,
    subjects_dir=subjects_dir,
    meg="helmet",
    surfaces=["white", "head"],
)