pooja debug.py

import os.path as op
import matplotlib.pyplot as plt
import mne
from mne.datasets import sample
from mne import setup_volume_source_space
from mne import make_forward_solution
from mne.minimum_norm import make_inverse_operator, apply_inverse
from nilearn.plotting import plot_stat_map
from nilearn.image import index_img

from joblib import Memory

mem = Memory('/tmp')

# Set dir
data_path = sample.data_path()
subject = 'sample'
data_dir = op.join(data_path, 'MEG', subject)
subjects_dir = op.join(data_path, 'subjects')
bem_dir = op.join(subjects_dir, subject, 'bem')
# Set file names
fname_aseg = op.join(subjects_dir, subject, 'mri', 'aseg.mgz')
fname_model = op.join(bem_dir, '%s-5120-bem.fif' % subject)
fname_bem = op.join(bem_dir, '%s-5120-bem-sol.fif' % subject)
fname_evoked = data_dir + '/sample_audvis-ave.fif'
fname_trans = data_dir + '/sample_audvis_raw-trans.fif'
fname_fwd = data_dir + '/sample_audvis-meg-oct-6-mixed-fwd.fif'
fname_cov = data_dir + '/sample_audvis-shrunk-cov.fif'
#########################################################################
labels_vol = ['Left-Hippocampus']  # labels_vol = ['Left-Hippocampus', ‘Right-Hippocampus’]
# labels_vol = ['Left-Cerebellum-Cortex', 'Right-Cerebellum-Cortex']
labels_vol = ['Left-Hippocampus', 'Right-Hippocampus']
# labels_vol = ['Left-Hippocampus', 'Right-Hippocampus',
#               'Left-Cerebellum-Cortex', 'Right-Cerebellum-Cortex']
sphere = (0, 0, 0, 120)
src = mem.cache(setup_volume_source_space)(subject, mri=fname_aseg, pos=2.0, sphere=sphere,
                                        volume_label=labels_vol, subjects_dir=subjects_dir, 
                                        verbose=True)
fwd = mem.cache(make_forward_solution)(fname_evoked, fname_trans, src, fname_bem, mindist=5.0,  
                            meg=True, eeg=False, n_jobs=1)
src_fwd = fwd['src']
n = sum(src_fwd[i]['nuse'] for i in range(len(src_fwd)))
print('the fwd src space contains %d spaces and %d points' % (len(src_fwd), n))
# Compute inverse solution 
snr = 3.0           
inv_method = 'MNE' # or MNE
lambda2 = 1.0 / snr ** 2
# Load data
condition = 'Left Auditory'
evoked = mne.read_evokeds(fname_evoked, condition=condition,
                          baseline=(None, 0))
noise_cov = mne.read_cov(fname_cov)
# Compute inverse operator
inverse_operator = make_inverse_operator(evoked.info, fwd, noise_cov,
                                         depth=None, fixed=False)
stcs = apply_inverse(evoked, inverse_operator, lambda2, inv_method,
                     pick_ori=None)
######################################################################
src = fwd['src']
t_ = stcs.times[200]

img = stcs.as_volume(src, mri_resolution=False)  # set True for full MRI resolution
t1_fname = data_path + '/subjects/sample/mri/T1.mgz'
# Plotting with nilearn ######################################################
threshold = 2*1e-9
# threshold = 1*1e-12
plot_stat_map(index_img(img, 200), t1_fname, threshold=threshold,
              title='%s (t=%.1f s.)' % (inv_method, t_),
              cut_coords=(-8, -35, -9))
plt.show()