plot_compare_morphing_options.py

# Authors: Denis A. Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)
import os
import os.path as op
import shlex
from subprocess import call

import numpy as np
import matplotlib.pyplot as plt

from mne import io
import mne

from mne.datasets import sample
from mne.minimum_norm import make_inverse_operator, apply_inverse
print(__doc__)

"""paths and params
"""

data_path = sample.data_path()

subject_from = 'sample'
subject_to = 'fsaverage'
subjects_dir = data_path + '/subjects'
os.environ.update(SUBJECTS_DIR=subjects_dir)

raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
event_id, tmin, tmax = 1, -0.2, 0.5

pick_ori = 'normal'
# pick_ori = None
ico = 5

"""Prepare data
"""
# Setup for reading the raw data
raw = io.Raw(raw_fname)
events = mne.read_events(event_fname)

# Set up pick list: EEG + MEG - bad channels (modify to your needs)
raw.info['bads'] += ['MEG 2443', 'EEG 053']  # bads + 2 more
picks = mne.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,
                       exclude='bads')

# Read epochs
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=(None, 0), preload=True,
                    reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))

evoked = epochs.average()  # average epochs to get the evoked response

"""pre-morphing (morphing the source space)
"""
call(shlex.split(
    'mne_setup_source_space --morph sample --ico %i '
    '--subject fsaverage --overwrite' % ico))

src_pre_morphed_name = op.join(
    subjects_dir, subject_from, 'bem', 'sample-fsaverage-ico-%s-src.fif' % ico)

src_morphed = mne.read_source_spaces(src_pre_morphed_name)

bem = mne.read_bem_solution(
    op.join(subjects_dir, subject_from, 'bem', 'sample-5120-bem-sol.fif'))

fwd = mne.make_forward_solution(
    evoked.info, 
    trans=op.join(
        data_path, 'MEG', subject_from, 'sample_audvis_raw-trans.fif'),
    src=src_morphed,
    n_jobs=8,
    bem=bem, fname=op.join(data_path , 'MEG', subject_from,
                           '%s-ico5-morphed-fwd.fif' % subject_from),
    overwrite=True)

fwd = mne.convert_forward_solution(fwd, surf_ori=True, force_fixed=False)

noise_cov = mne.read_cov(
    op.join(data_path, 'MEG', subject_from, 'sample_audvis-shrunk-cov.fif'))

invop = make_inverse_operator(
    evoked.info, fwd, noise_cov, loose=0.2, depth=0.8, fixed=False)

snr = 3.0
lambda2 = 1.0 / snr ** 2

stc_pre_morphed = apply_inverse(
    evoked, invop, lambda2, "dSPM", pick_ori=pick_ori)


"""other options: morphing the stc"""

fname_fwd_meg = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
fwd_from = mne.read_forward_solution(fname_fwd_meg)
fwd_from = mne.convert_forward_solution(
    # fwd_from, surf_ori=True, force_fixed=True)
    fwd_from, surf_ori=True, force_fixed=False)

invop_from = make_inverse_operator(
    evoked.info, fwd_from, noise_cov, loose=0.2, depth=0.8, fixed=False)

stc_from = apply_inverse(
    evoked, invop_from, lambda2, "dSPM", pick_ori=pick_ori)
vertices_to = [np.arange(10242), np.arange(10242)]
stc_to = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=1,
                        grade=vertices_to, subjects_dir=subjects_dir)
stc_to.save('%s_audvis-meg' % subject_to)

morph_mat = mne.compute_morph_matrix(subject_from, subject_to,
                                     stc_from.vertices, vertices_to,
                                     subjects_dir=subjects_dir)
stc_to_2 = mne.morph_data_precomputed(subject_from, subject_to,
                                      stc_from, vertices_to, morph_mat)
stc_to_2.save('%s_audvis-meg_2' % subject_to)


fig = plt.figure()
plt.plot(stc_from.times, stc_from.data.mean(axis=0), 'r', label='from')
plt.plot(stc_to.times, stc_to.data.mean(axis=0), 'b', label='to')
plt.plot(stc_to_2.times, stc_to_2.data.mean(axis=0), 'g', label='to_2')
plt.plot(
    stc_to_2.times, stc_pre_morphed.data.mean(axis=0), 'y',
    label='pre_morphed')
plt.xlabel('time (ms)')
plt.ylabel('Mean Source amplitude')
plt.legend()
plt.show()
fig.savefig('morphing-options_pick-ori-%s_ico-%i.png' % (pick_ori, ico))

brain = stc_pre_morphed.plot(hemi='both', clim='auto')
brain.set_time(100)
brain.save_image('stc_pre-morphed_pick-ori-%s_ico-%i.png' % (pick_ori, ico))
brain.close()
brain = stc_to.plot(hemi='both', clim='auto')
brain.set_time(100)
brain.save_image('stc_post-morphed_pick-ori-%s_ico-%i.png' % (pick_ori, ico))
brain.close()

"""in label"""

aud_lh = mne.read_label(data_path + '/MEG/sample' + '/labels/Aud-lh.label')
aud_rh = mne.read_label(data_path + '/MEG/sample' + '/labels/Aud-rh.label')


fig, axes = plt.subplots(1, 2, figsize=(12, 6), sharex=True, sharey=True)
ax1, ax2 = axes.ravel()
plt.suptitle('In auditory labels')
times = stc_to.times * 1e3
ax1.plot(
    times, stc_from.in_label(aud_rh).data.mean(axis=0), color='steelblue',
    label='native')
ax1.plot(
    times, stc_pre_morphed.in_label(aud_rh).data.mean(axis=0), color='red',
    label='morphed<-fsaverage')
ax2.plot(times, stc_from.in_label(aud_lh).data.mean(axis=0), color='steelblue',
    label='native')
ax2.plot(
    times, stc_pre_morphed.in_label(aud_lh).data.mean(axis=0), color='red',
    label='morphed<-fsaverage')

ax1.set_title('auditory rh')
ax2.set_title('auditory lh')
ax1.set_ylabel('dSPM source estimate [Arn. U.]')
ax2.set_xlabel('Time [ms]')
plt.legend(loc='lower right')
fig.savefig('morphing-in-label_pick-ori-%s_ico-%i.png' % (pick_ori, ico))