sample_mri_pos.py

sample_mri_pos.py

# Create a DigMontage for sample in MRI coordinates

import os.path as op
import numpy as np
import nibabel
import mne
from mne.transforms import apply_trans
from mne.io.constants import FIFF
from mne.channels import make_dig_montage, read_custom_montage

# Create:
#
# - ``eeg_vox``: EEG electrodes in RAS voxel coordinates
# - ``fid_vox``: fiducial (LPA, nasion, RPA) locations in RAS voxel coordinates

data_path = mne.datasets.sample.data_path()
subjects_dir = op.join(data_path, 'subjects')
fname_T1 = op.join(subjects_dir, 'sample', 'mri', 'T1.mgz')
trans = mne.read_trans(
    op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw-trans.fif'))
raw = mne.io.read_raw_fif(
    op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif'))

# Get fiducials and EEG locations in MRI voxel coords
trans_mm = trans['trans'].copy()
trans_mm[:3, 3] *= 1e3
img = nibabel.load(fname_T1)
vox2mri = img.header.get_vox2ras_tkr()
head_to_vox = np.dot(np.linalg.inv(vox2mri), trans_mm)
fid_vox = apply_trans(head_to_vox, 1000 * np.array([d['r'] for d in sorted(
    [d for d in raw.info['dig'] if d['kind'] == FIFF.FIFFV_POINT_CARDINAL],
    key=lambda d: d['ident'])]))
eeg_vox = apply_trans(head_to_vox, 1000 * np.array(
    [d['r'] for d in raw.info['dig'] if d['kind'] == FIFF.FIFFV_POINT_EEG]))
img, affine, header = img.get_data().astype(float), img.affine, img.header
# Mark them in the MRI with 5x5x5mm cubes
for coord in np.round(np.concatenate([fid_vox, eeg_vox])).astype(int):
    img[tuple(slice(c - 2, c + 2) for c in coord)] = 250
img = nibabel.freesurfer.mghformat.MGHImage(img, affine, header=header)
nibabel.save(img, 'T1_electrodes.mgz')
# change "raw" to be more like what you'd acquire: no dig, no channel locations
raw.load_data().pick_types(meg=False, eeg=True, exclude=())
del head_to_vox, trans

# Get our MRI channel locations first into RAS space, in meters.
img = nibabel.load(fname_T1)
vox2mri = img.header.get_vox2ras_tkr()  # MRI voxel -> FreeSurfer RAS (mm)
vox2mri[:3] /= 1000.  # mm -> m
fid_mri = apply_trans(vox2mri, fid_vox)
eeg_mri = apply_trans(vox2mri, eeg_vox)

# make our digitization object for sanity checking

eeg_names = ['EEG 000'] + [raw.ch_names[pick] for pick in mne.pick_types(
    raw.info, meg=False, eeg=True, exclude=())]
assert len(eeg_names) == len(eeg_mri)
ch_pos = {ch_name: pos for ch_name, pos in zip(eeg_names, eeg_mri)}
lpa, nasion, rpa = fid_mri
dig_montage = make_dig_montage(ch_pos, nasion, lpa, rpa, coord_frame='mri')
dig_montage.plot()

# Output an ASA electrode file
fname = 'sample_mri_montage.elc'
all_names = ['LPA', 'Nz', 'RPA'] + eeg_names
all_pos = np.concatenate([fid_mri, eeg_mri], axis=0)
assert len(all_names) == len(all_pos)
with open(fname, 'w') as fid:
    fid.write('ReferenceLabel\navg\nUnitPosition\tmm\n'
              'NumberPositions=\t%s\nPositions\n' % (len(all_names),))
    for p in all_pos:
        fid.write('%0.4f %0.4f %0.4f\n' % tuple(1000 * p))
    fid.write('Labels\n')
    for name in all_names:
        fid.write('%s\n' % (name,))
dig_montage2 = read_custom_montage(fname, head_size=None, coord_frame='mri')
dig_montage2.plot()