plot_topo_grad_complex.py

# Author: Jean-Remi King <[email protected]>
#
# License: BSD (3-clause)

import numpy as np
import matplotlib.pyplot as plt

from mne import read_evokeds
from mne.datasets.megsim import load_data
from mne.viz.topomap import _griddata, _prepare_topo_plot

# Preprocess data
epochs_fnames = load_data(condition='visual', data_format='single-trial',
                          data_type='simulation')

# Take only 10 trials from the same simulation setup.
epochs_fnames = [f for f in epochs_fnames if 'sim6_trial_' in f][:10]
evokeds = [read_evokeds(f)[0] for f in epochs_fnames]
evoked = sum(evokeds[1:], evokeds[0])

# Select time and chan
data = evoked.data[:, np.where(evoked.times >= .350)[0][0]]

# Construct topo
picks, pos, _, _, ch_type = _prepare_topo_plot(evoked, 'grad', None)
pos_x, pos_y = pos[:, 0], pos[:, 1]
xmin, xmax = pos_x.min(), pos_x.max()
ymin, ymax = pos_y.min(), pos_y.max()

res = 64
xi = np.linspace(xmin, xmax, res)
yi = np.linspace(ymin, ymax, res)
Xi, Yi = np.meshgrid(xi, yi)
Zi1 = _griddata(pos_x, pos_y, data[picks[::2]], Xi, Yi)
Zi2 = _griddata(pos_x, pos_y, data[picks[1::2]], Xi, Yi)


def rgbscale(x, m=-3 * 1e-12, M=3 * 1e-12):
    [ni, nj, nk] = x.shape
    for i in range(ni):
        for j in range(nj):
            for k in range(nk):
                if x[i, j, k] < m:
                    x[i, j, k] = 0.
                elif x[i, j, k] > M:
                    x[i, j, k] = 1.
                else:
                    x[i, j, k] += m
                    x[i, j, k] /= (M - m)
    return x


Zi = np.concatenate((Zi1[..., None],
                     np.zeros(Zi1.shape)[..., None],
                     Zi2[..., None]), axis=2)
Zi = rgbscale(Zi)

im = plt.imshow(Zi, origin='lower', aspect='equal', interpolation='bilinear',
                extent=(xmin, xmax, ymin, ymax))