debug_src_vol_connectivity.py

import numpy as np
from scipy import sparse
import matplotlib.pyplot as plt

import mne
from mne.datasets import sample
from mne import read_evokeds
from mne.minimum_norm import read_inverse_operator

print(__doc__)

data_path = sample.data_path()
fname_inv = data_path + '/MEG/sample/sample_audvis-meg-vol-7-meg-inv.fif'
fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif'
subject = 'sample'

snr = 3.0
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

# Load data
evoked = read_evokeds(fname_evoked, condition=0, baseline=(None, 0))
inverse_operator = read_inverse_operator(fname_inv)
src = inverse_operator['src']


def spatial_src_connectivity(src, subject=None):
    # Export result as a 4D nifti object
    from sklearn.feature_extraction import grid_to_graph

    assert src[0]['type'] == 'vol'

    if subject is None:
        subject = src[0]['subject_his_id']

    vertices = np.where(src[0]['inuse'])[0]
    n_vertices = len(vertices)
    data = (1 + np.arange(n_vertices))[:, np.newaxis]
    stc_tmp = mne.VolSourceEstimate(data, vertices, tmin=0., tstep=1.,
                                    subject=subject)
    img = stc_tmp.as_volume(src, mri_resolution=False)
    img_data = img.get_data()[:, :, :, 0]
    img_data = img_data.T
    graph = grid_to_graph(*img_data.shape, mask=(img_data != 0))
    return graph

graph = spatial_src_connectivity(src)
print(graph.shape)


stc_data = np.zeros((graph.shape[0], 1))
stc_data[2000] = 1.


def mesh_smoothing_matrix(A):
    inv_conn_idx = .75 / A.sum(axis=1)
    A = sparse.spdiags(inv_conn_idx.T, 0, *A.shape).dot(A)
    A.setdiag(np.ones(A.shape[0]))
    return A

A = mesh_smoothing_matrix(graph)

for i in range(15):
    stc_data = A.dot(stc_data)

vertices = np.where(src[0]['inuse'])[0]
stc_tmp = mne.VolSourceEstimate(stc_data, vertices, tmin=0., tstep=1.,
                                subject=subject)
img = stc_tmp.as_volume(src, mri_resolution=False)

img_data = img.get_data()[:, :, :, 0]

plt.imshow(img_data[15])
plt.show()