From adce5881b814bef12e0686a801eaadd81b29d507 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Mon, 11 May 2020 21:01:33 -0400 Subject: [PATCH 01/27] Adding updated tutorial with ecog data. --- mne/datasets/utils.py | 5 +++-- tutorials/misc/plot_ecog.py | 27 ++++++++++++++++++++++++++- 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 7f9131b4fc7..38c353e9b40 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -254,8 +254,9 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, bst_resting='https://osf.io/m7bd3/download?version=3'), fake='https://github.com/mne-tools/mne-testing-data/raw/master/' 'datasets/foo.tgz', - misc='https://codeload.github.com/mne-tools/mne-misc-data/' - 'tar.gz/%s' % releases['misc'], + # misc='https://codeload.github.com/mne-tools/mne-misc-data/' + # 'tar.gz/%s' % releases['misc'], + misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', sample='https://osf.io/86qa2/download?version=5', somato='https://osf.io/tp4sg/download?version=6', spm='https://osf.io/je4s8/download?version=2', diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 1d1b57d6446..111ee45010b 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -9,6 +9,7 @@ """ # Authors: Eric Larson # Chris Holdgraf +# Edited: Adam Li # # License: BSD (3-clause) @@ -28,7 +29,8 @@ mat = loadmat(mne.datasets.misc.data_path() + '/ecog/sample_ecog.mat') ch_names = mat['ch_names'].tolist() elec = mat['elec'] # electrode positions given in meters -# Now we make a montage stating that the sEEG contacts are in head + +# Now we make a montage stating that the ECoG contacts are in head # coordinate system (although they are in MRI). This is compensated # by the fact that below we do not specicty a trans file so the Head<->MRI # transform is the identity. @@ -77,3 +79,26 @@ ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') ax.set_axis_off() plt.show() + +############################################################################### +# Sometimes it is useful to create an animation of the ECoG activity over time. +# Wee can visualize + +# We'll once again plot the surface, then take a snapshot. +fig_scatter = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, + surfaces='pial') +mne.viz.set_3d_view(fig_scatter, 200, 70) +xy, im = snapshot_brain_montage(fig_scatter, montage) + +# Convert from a dictionary to array to plot +xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) + +# Define an arbitrary "activity" pattern for viz +activity = np.linspace(100, 200, xy_pts.shape[0]) + +# This allows us to use matplotlib to create arbitrary 2d scatterplots +_, ax = plt.subplots(figsize=(10, 10)) +ax.imshow(im) +ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') +ax.set_axis_off() +plt.show() From 502307e915e4df1192f7775c87cca70e8ac927de Mon Sep 17 00:00:00 2001 From: Adam Li Date: Mon, 11 May 2020 23:52:14 -0400 Subject: [PATCH 02/27] Trying to update tutorial w/ misc data. --- mne/datasets/utils.py | 6 ++-- tutorials/misc/plot_ecog.py | 64 ++++++++++++++++++++++++++----------- 2 files changed, 50 insertions(+), 20 deletions(-) diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 38c353e9b40..ac2ecf8768d 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -278,7 +278,8 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, archive_names = dict( fieldtrip_cmc='SubjectCMC.zip', kiloword='MNE-kiloword-data.tar.gz', - misc='mne-misc-data-%s.tar.gz' % releases['misc'], + # misc='mne-misc-data-%s.tar.gz' % releases['misc'], + misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', mtrf='mTRF_1.5.zip', multimodal='MNE-multimodal-data.tar.gz', fnirs_motor='MNE-fNIRS-motor-data.tgz', @@ -294,7 +295,8 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, # original folder names that get extracted (only needed if the # archive does not extract the right folder name; e.g., usually GitHub) folder_origs = dict( # not listed means None (no need to move) - misc='mne-misc-data-%s' % releases['misc'], + # misc='mne-misc-data-%s' % releases['misc'], + misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', testing='mne-testing-data-%s' % releases['testing'], ) # finally, where we want them to extract to (only needed if the folder name diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 111ee45010b..886cf040818 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -73,32 +73,60 @@ # Define an arbitrary "activity" pattern for viz activity = np.linspace(100, 200, xy_pts.shape[0]) -# This allows us to use matplotlib to create arbitrary 2d scatterplots -_, ax = plt.subplots(figsize=(10, 10)) -ax.imshow(im) -ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') -ax.set_axis_off() -plt.show() +# # This allows us to use matplotlib to create arbitrary 2d scatterplots +# _, ax = plt.subplots(figsize=(10, 10)) +# ax.imshow(im) +# ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') +# ax.set_axis_off() +# plt.show() ############################################################################### # Sometimes it is useful to create an animation of the ECoG activity over time. -# Wee can visualize +# We can visualize say the gamma frequency of the ECoG activity on the brain +# using MNE functions. -# We'll once again plot the surface, then take a snapshot. -fig_scatter = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, - surfaces='pial') -mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, montage) +# first we'll load in the sample dataset +raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf') -# Convert from a dictionary to array to plot -xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) +# attach montage +raw.set_montage(montage, on_missing='warn') + +# perform gamma band frequency +tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(raw, freqs=np.linspace(30, 90, 60), + n_cycles=7) # Define an arbitrary "activity" pattern for viz -activity = np.linspace(100, 200, xy_pts.shape[0]) +activity = tfr_pwr.data.mean(axis=1) -# This allows us to use matplotlib to create arbitrary 2d scatterplots -_, ax = plt.subplots(figsize=(10, 10)) +# create animation over the entire time period of 10 seconds +# from celluloid import Camera +import matplotlib.animation as animation +fig, ax = plt.subplots(figsize=(10, 10)) +# camera = Camera(fig) ax.imshow(im) -ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') ax.set_axis_off() + +paths = ax.scatter([], c=[], s=200, cmap='coolwarm') + +# initialization function +def init(): + # creating an empty plot/frame + paths.set_data([], []) + return paths, + +# animation function +def animate(i): + paths = ax.scatter(*xy_pts.T, c=activity[:, i], s=200, cmap='coolwarm') + + # appending new points to x, y axes points list + # line.set_data(xdata, ydata) + return paths, + +# call the animator +anim = animation.FuncAnimation(fig, animate, + init_func=init, + frames=500, interval=20, blit=True) + + +animation = camera.animate() plt.show() From 5c1e794dc06a3d129f346bc4968102af25cfd9e9 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 00:56:48 -0400 Subject: [PATCH 03/27] Adding prelim ecog plotting. --- tutorials/misc/plot_ecog.py | 96 +++++++++++++++++++++++-------------- 1 file changed, 60 insertions(+), 36 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 886cf040818..d51b7e789ed 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -14,6 +14,12 @@ # License: BSD (3-clause) import numpy as np +from sys import platform as sys_pf +# if sys_pf == 'darwin': +# import matplotlib +# matplotlib.use("Qt5Agg") +import matplotlib +matplotlib.use("macOSX") import matplotlib.pyplot as plt from scipy.io import loadmat @@ -30,12 +36,21 @@ ch_names = mat['ch_names'].tolist() elec = mat['elec'] # electrode positions given in meters +from mne_bids.tsv_handler import _from_tsv +elec_tsv = _from_tsv(mne.datasets.misc.data_path() + '/ecog/sample_ecog_electrodes.tsv') +ch_names = elec_tsv['name'] +ch_coords = np.vstack((elec_tsv['x'], elec_tsv['y'], elec_tsv['z'])).T.astype(float) +ch_pos = dict(zip(ch_names, ch_coords)) +montage = mne.channels.make_dig_montage(ch_pos, + coord_frame='head') +print(ch_names) +print(ch_pos) # Now we make a montage stating that the ECoG contacts are in head # coordinate system (although they are in MRI). This is compensated # by the fact that below we do not specicty a trans file so the Head<->MRI # transform is the identity. -montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)), - coord_frame='head') +# montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)), +# coord_frame='head') print('Created %s channel positions' % len(ch_names)) ############################################################################### @@ -74,10 +89,10 @@ activity = np.linspace(100, 200, xy_pts.shape[0]) # # This allows us to use matplotlib to create arbitrary 2d scatterplots -# _, ax = plt.subplots(figsize=(10, 10)) -# ax.imshow(im) -# ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') -# ax.set_axis_off() +_, ax = plt.subplots(figsize=(10, 10)) +ax.imshow(im) +ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') +ax.set_axis_off() # plt.show() ############################################################################### @@ -88,45 +103,54 @@ # first we'll load in the sample dataset raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf') +# drop bad channels +raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) + # attach montage raw.set_montage(montage, on_missing='warn') # perform gamma band frequency -tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(raw, freqs=np.linspace(30, 90, 60), - n_cycles=7) - +epoch = mne.EpochsArray(raw.get_data()[np.newaxis, ...], info=raw.info) +print(epoch) +# print(epoch.shape) +tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(30, 90, 60), + n_cycles=3) +print(tfr_pwr) # Define an arbitrary "activity" pattern for viz -activity = tfr_pwr.data.mean(axis=1) +gamma_activity = tfr_pwr.data.mean(axis=(1, 2)) -# create animation over the entire time period of 10 seconds -# from celluloid import Camera -import matplotlib.animation as animation -fig, ax = plt.subplots(figsize=(10, 10)) -# camera = Camera(fig) -ax.imshow(im) -ax.set_axis_off() - -paths = ax.scatter([], c=[], s=200, cmap='coolwarm') - -# initialization function -def init(): - # creating an empty plot/frame - paths.set_data([], []) - return paths, +tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(1, 30, 60), + n_cycles=3) +low_activity = tfr_pwr.data.mean(axis=(1, 2)) -# animation function -def animate(i): - paths = ax.scatter(*xy_pts.T, c=activity[:, i], s=200, cmap='coolwarm') - # appending new points to x, y axes points list - # line.set_data(xdata, ydata) - return paths, +_, ax = plt.subplots(figsize=(10, 10)) +# show activity between low frequency and higher frequencies +# We'll once again plot the surface, then take a snapshot. +fig_scatter = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, + surfaces='pial') +mne.viz.set_3d_view(fig_scatter, 200, 70) +xy, im = snapshot_brain_montage(fig_scatter, montage) +# Convert from a dictionary to array to plot +xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) -# call the animator -anim = animation.FuncAnimation(fig, animate, - init_func=init, - frames=500, interval=20, blit=True) +ax.imshow(im) +ax.set_axis_off() +ax.scatter(*xy_pts.T, c=gamma_activity, s=200, cmap='coolwarm') +# ax.set_title("Gamma frequency (30-90 Hz)") +_, ax = plt.subplots(figsize=(10, 10)) +# show activity between low frequency and higher frequencies +# We'll once again plot the surface, then take a snapshot. +fig_scatter = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, + surfaces='pial') +mne.viz.set_3d_view(fig_scatter, 200, 70) +xy, im = snapshot_brain_montage(fig_scatter, montage) +# Convert from a dictionary to array to plot +xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) +ax.imshow(im) +ax.set_axis_off() +ax.scatter(*xy_pts.T, c=low_activity, s=200, cmap='coolwarm') +# ax.set_title("Low frequency (0-30 Hz)") -animation = camera.animate() plt.show() From 068e48713f2e3ea9e06b973f0c2fb0b06cb153b9 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Tue, 12 May 2020 09:56:26 -0400 Subject: [PATCH 04/27] FIX: Fetching --- mne/datasets/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index ac2ecf8768d..d27deb21332 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -279,7 +279,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, fieldtrip_cmc='SubjectCMC.zip', kiloword='MNE-kiloword-data.tar.gz', # misc='mne-misc-data-%s.tar.gz' % releases['misc'], - misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', + misc='mne-misc-data.tar.gz', mtrf='mTRF_1.5.zip', multimodal='MNE-multimodal-data.tar.gz', fnirs_motor='MNE-fNIRS-motor-data.tgz', @@ -296,7 +296,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, # archive does not extract the right folder name; e.g., usually GitHub) folder_origs = dict( # not listed means None (no need to move) # misc='mne-misc-data-%s' % releases['misc'], - misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', + misc='adam2392-mne-misc-data-75968e3', testing='mne-testing-data-%s' % releases['testing'], ) # finally, where we want them to extract to (only needed if the folder name @@ -321,7 +321,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, bst_raw='fa2efaaec3f3d462b319bc24898f440c', bst_resting='70fc7bf9c3b97c4f2eab6260ee4a0430'), fake='3194e9f7b46039bb050a74f3e1ae9908', - misc='84e606998ac379ef53029b3b1cf37918', + misc='a82a6888d4266d3c6d5aacc8f305b5b8', sample='12b75d1cb7df9dfb4ad73ed82f61094f', somato='ea825966c0a1e9b2f84e3826c5500161', spm='9f43f67150e3b694b523a21eb929ea75', From 7fdb1bc3d6e3a6f222eee9d881a9beafb0a04104 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Tue, 12 May 2020 10:08:45 -0400 Subject: [PATCH 05/27] FIX: Fix example --- tutorials/misc/plot_ecog.py | 46 +++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 886cf040818..8f8674ee5e2 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -15,6 +15,7 @@ import numpy as np import matplotlib.pyplot as plt +import matplotlib.animation as animation from scipy.io import loadmat import mne @@ -22,11 +23,14 @@ print(__doc__) +misc_path = mne.datasets.misc.data_path() +sample_path = mne.datasets.sample.data_path() + ############################################################################### # Let's load some ECoG electrode locations and names, and turn them into # a :class:`mne.channels.DigMontage` class. -mat = loadmat(mne.datasets.misc.data_path() + '/ecog/sample_ecog.mat') +mat = loadmat(misc_path + '/ecog/sample_ecog.mat') ch_names = mat['ch_names'].tolist() elec = mat['elec'] # electrode positions given in meters @@ -50,7 +54,7 @@ # .. note:: These are not real electrodes for this subject, so they # do not align to the cortical surface perfectly. -subjects_dir = mne.datasets.sample.data_path() + '/subjects' +subjects_dir = sample_path + '/subjects' fig = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, surfaces=['pial']) mne.viz.set_3d_view(fig, 200, 70) @@ -86,47 +90,45 @@ # using MNE functions. # first we'll load in the sample dataset -raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf') +raw = mne.io.read_raw_edf(misc_path + '/ecog/sample_ecog.edf') # attach montage -raw.set_montage(montage, on_missing='warn') +raw.set_montage(montage, on_missing='ignore') # perform gamma band frequency -tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(raw, freqs=np.linspace(30, 90, 60), - n_cycles=7) +events = [[raw.first_samp, 0, 1]] +epochs = mne.Epochs(raw, events, tmin=0, tmax=raw.times[-1], baseline=None, + preload=True) +tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet( + epochs, freqs=np.linspace(30, 90, 60), n_cycles=7) # Define an arbitrary "activity" pattern for viz activity = tfr_pwr.data.mean(axis=1) # create animation over the entire time period of 10 seconds -# from celluloid import Camera -import matplotlib.animation as animation fig, ax = plt.subplots(figsize=(10, 10)) -# camera = Camera(fig) ax.imshow(im) ax.set_axis_off() +vmin, vmax = np.percentile(activity, [10, 90]) + +paths = ax.scatter(*xy_pts.T, c=np.zeros(len(xy_pts)), s=200, cmap='plasma_r', + vmin=vmin, vmax=vmax) -paths = ax.scatter([], c=[], s=200, cmap='coolwarm') # initialization function def init(): - # creating an empty plot/frame - paths.set_data([], []) + """Create an empty frame.""" return paths, + # animation function def animate(i): - paths = ax.scatter(*xy_pts.T, c=activity[:, i], s=200, cmap='coolwarm') - - # appending new points to x, y axes points list - # line.set_data(xdata, ydata) + """Animate the plot.""" + paths.set_array(activity[:, i]) return paths, -# call the animator -anim = animation.FuncAnimation(fig, animate, - init_func=init, - frames=500, interval=20, blit=True) - -animation = camera.animate() +# call the animator +anim = animation.FuncAnimation( + fig, animate, init_func=init, frames=500, interval=20, blit=True) plt.show() From 5e42269ca78fbc142491df8c53ec7e078ffbe298 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 10:42:09 -0400 Subject: [PATCH 06/27] Merging. --- Untitled.ipynb | 277 ++++++++++++++++++++++++++++++++++++ tutorials/misc/plot_ecog.py | 138 +++++++++--------- 2 files changed, 348 insertions(+), 67 deletions(-) create mode 100644 Untitled.ipynb diff --git a/Untitled.ipynb b/Untitled.ipynb new file mode 100644 index 00000000000..667123ca79c --- /dev/null +++ b/Untitled.ipynb @@ -0,0 +1,277 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib\n", + "# matplotlib.use(\"macOSX\")\n", + "import matplotlib.pyplot as plt\n", + "from scipy.io import loadmat\n", + "\n", + "from sklearn.preprocessing import StandardScaler\n", + "\n", + "import mne\n", + "from mne.viz import plot_alignment, snapshot_brain_montage" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9', 'C10', 'C11', 'C12', 'C13', 'C14', 'C15', 'C16', 'C17', 'C18', 'C19', 'C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C27', 'C28', 'C29', 'C30', 'C31', 'C32', 'C33', 'C34', 'C35', 'C36', 'C37', 'C38', 'C39', 'C40', 'C41', 'C42', 'C43', 'C44', 'C45', 'C46', 'C47', 'C48', 'C49', 'C50', 'C51', 'C52', 'C53', 'C54', 'C55', 'C56', 'C57', 'C58', 'C59', 'C60', 'C61', 'C62', 'C63', 'C64']\n", + "{'C1': array([-0.0132, -0.0654, 0.0702]), 'C2': array([-0.0145, -0.0594, 0.0781]), 'C3': array([-0.0156, -0.0518, 0.0838]), 'C4': array([-0.0154, -0.0453, 0.0895]), 'C5': array([-0.0149, -0.0385, 0.094 ]), 'C6': array([-0.0131, -0.029 , 0.0985]), 'C7': array([-0.0136, -0.0203, 0.1009]), 'C8': array([-0.0109, -0.0093, 0.1035]), 'C9': array([-0.0215, -0.0664, 0.0691]), 'C10': array([-0.0233, -0.0607, 0.0769]), 'C11': array([-0.0247, -0.0534, 0.0822]), 'C12': array([-0.0263, -0.0461, 0.0871]), 'C13': array([-0.0261, -0.0383, 0.0918]), 'C14': array([-0.0253, -0.0299, 0.0949]), 'C15': array([-0.0233, -0.0193, 0.0989]), 'C16': array([-0.0211, -0.008 , 0.1026]), 'C17': array([-0.0304, -0.0674, 0.0678]), 'C18': array([-0.0328, -0.061 , 0.0738]), 'C19': array([-0.0337, -0.0539, 0.0793]), 'C20': array([-0.0344, -0.0466, 0.0842]), 'C21': array([-0.0363, -0.0388, 0.0878]), 'C22': array([-0.0345, -0.0288, 0.0926]), 'C23': array([-0.0335, -0.018 , 0.0955]), 'C24': array([-0.0299, -0.008 , 0.0977]), 'C25': array([-0.0395, -0.0661, 0.0632]), 'C26': array([-0.0417, -0.0599, 0.069 ]), 'C27': array([-0.0426, -0.0529, 0.0744]), 'C28': array([-0.0438, -0.0456, 0.0795]), 'C29': array([-0.0457, -0.0376, 0.083 ]), 'C30': array([-0.0446, -0.0273, 0.0881]), 'C31': array([-0.043 , -0.0176, 0.0901]), 'C32': array([-0.0396, -0.0063, 0.0884]), 'C33': array([-0.0446, -0.0657, 0.0581]), 'C34': array([-0.0482, -0.0593, 0.0632]), 'C35': array([-0.0505, -0.0523, 0.0694]), 'C36': array([-0.0536, -0.0438, 0.0731]), 'C37': array([-0.0548, -0.0357, 0.0769]), 'C38': array([-0.0533, -0.0254, 0.0808]), 'C39': array([-0.0512, -0.0156, 0.0823]), 'C40': array([-0.0473, -0.0044, 0.0808]), 'C41': array([-0.0518, -0.0637, 0.0513]), 'C42': array([-0.055 , -0.0572, 0.0574]), 'C43': array([-0.0575, -0.0501, 0.0625]), 'C44': array([-0.0606, -0.0414, 0.0666]), 'C45': array([-0.0614, -0.0334, 0.0709]), 'C46': array([-0.0612, -0.023 , 0.0733]), 'C47': array([-0.0589, -0.0135, 0.0743]), 'C48': array([-0.0553, -0.0025, 0.0739]), 'C49': array([-0.0575, -0.0605, 0.0455]), 'C50': array([-0.0615, -0.0536, 0.05 ]), 'C51': array([-0.0635, -0.0474, 0.0554]), 'C52': array([-0.066 , -0.0385, 0.0584]), 'C53': array([-0.0671, -0.0306, 0.063 ]), 'C54': array([-0.0672, -0.0204, 0.0644]), 'C55': array([-0.0643, -0.0112, 0.0638]), 'C56': array([-0.0602, -0.0007, 0.064 ]), 'C57': array([-0.0614, -0.0574, 0.0372]), 'C58': array([-0.0659, -0.05 , 0.0412]), 'C59': array([-0.0681, -0.0434, 0.0463]), 'C60': array([-0.0694, -0.0357, 0.0503]), 'C61': array([-0.0704, -0.0277, 0.053 ]), 'C62': array([-0.0705, -0.0177, 0.0541]), 'C63': array([-0.068 , -0.0074, 0.053 ]), 'C64': array([-0.0642, 0.0016, 0.0539])}\n" + ] + } + ], + "source": [ + "\n", + "mat = loadmat(mne.datasets.misc.data_path() + '/ecog/sample_ecog.mat')\n", + "ch_names = mat['ch_names'].tolist()\n", + "elec = mat['elec'] # electrode positions given in meters\n", + "\n", + "from mne_bids.tsv_handler import _from_tsv\n", + "elec_tsv = _from_tsv(mne.datasets.misc.data_path() + '/ecog/sample_ecog_electrodes.tsv')\n", + "ch_names = elec_tsv['name']\n", + "ch_coords = np.vstack((elec_tsv['x'], elec_tsv['y'], elec_tsv['z'])).T.astype(float)\n", + "ch_pos = dict(zip(ch_names, ch_coords))\n", + "montage = mne.channels.make_dig_montage(ch_pos,\n", + " coord_frame='head')\n", + "print(ch_names)\n", + "print(ch_pos)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Extracting EDF parameters from /Users/adam2392/mne_data/MNE-misc-data/ecog/sample_ecog.edf...\n", + "EDF file detected\n", + "Setting channel info structure...\n", + "Creating raw.info structure...\n", + "1 matching events found\n", + "No baseline correction applied\n", + "Not setting metadata\n", + "0 projection items activated\n", + "0 bad epochs dropped\n", + "\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + ":19: RuntimeWarning: DigMontage is only a subset of info. There are 14 channel positions not present in the DigMontage. The required channels are: ['BTM1', 'BTM2', 'BTM3', 'BTM4', 'BTM5', 'BTM6', 'BTP1', 'BTP2', 'BTP3', 'BTP4', 'BTP5', 'BTP6', 'EKGL', 'EKGR']\n", + " raw.set_montage(montage, on_missing='warn')\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "(64,)\n" + ] + } + ], + "source": [ + "###############################################################################\n", + "# Now that we have our electrode positions in MRI coordinates, we can create\n", + "# our measurement info structure.\n", + "\n", + "info = mne.create_info(ch_names, 1000., 'ecog').set_montage(montage)\n", + "\n", + "###############################################################################\n", + "# Now that we have our electrode positions in MRI coordinates, we can load in\n", + "# our corresponding time-series data. We then compute a time-frequency\n", + "# representation of the data (i.e. 1-30, or 30-90 Hz).\n", + "\n", + "# first we'll load in the sample dataset\n", + "raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf')\n", + "\n", + "# drop bad channels\n", + "raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names])\n", + "\n", + "# attach montage\n", + "raw.set_montage(montage, on_missing='warn')\n", + "\n", + "# perform gamma band frequency\n", + "epoch = mne.EpochsArray(raw.get_data()[np.newaxis, ...], info=raw.info)\n", + "print(epoch)\n", + "# print(epoch.shape)\n", + "tfr_pwr, _ = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(30, 90, 60),\n", + " n_cycles=3)\n", + "print(tfr_pwr)\n", + "# Define an arbitrary \"activity\" pattern for viz\n", + "gamma_activity = tfr_pwr.data.mean(axis=(1, 2))\n", + "print(gamma_activity.shape)\n", + "gamma_activity = StandardScaler().fit_transform(gamma_activity[:, np.newaxis])\n", + "\n", + "tfr_pwr, _ = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(1, 30, 60),\n", + " n_cycles=3)\n", + "low_activity = tfr_pwr.data.mean(axis=(1, 2))\n", + "low_activity = StandardScaler().fit_transform(low_activity[:, np.newaxis])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0.30451935]\n", + " [ 0.35250178]\n", + " [ 0.18202274]\n", + " [ 0.4322955 ]\n", + " [ 0.40049513]\n", + " [-0.72390702]\n", + " [-0.72394051]\n", + " [-0.78456946]\n", + " [-0.36576045]\n", + " [-0.96004766]\n", + " [-0.97569081]\n", + " [-0.87451106]\n", + " [ 0.91036912]\n", + " [-0.31239582]\n", + " [-0.43099794]\n", + " [-0.5732088 ]\n", + " [ 0.53345374]\n", + " [-0.91170083]\n", + " [-0.32429372]\n", + " [-0.95974625]\n", + " [-0.56329589]\n", + " [ 2.94155747]\n", + " [-0.66383375]\n", + " [-0.50243725]\n", + " [ 3.00058893]\n", + " [-0.62034532]\n", + " [ 0.61977523]\n", + " [-0.2694325 ]\n", + " [ 2.1429474 ]\n", + " [-0.98527046]\n", + " [ 1.41722789]\n", + " [-0.29688988]\n", + " [-0.94528061]\n", + " [ 0.54137737]\n", + " [ 0.74136725]\n", + " [ 0.0561976 ]\n", + " [ 0.35904437]\n", + " [-0.76816171]\n", + " [-0.26112042]\n", + " [-0.56606875]\n", + " [ 1.82960786]\n", + " [-0.56582553]\n", + " [-0.96589049]\n", + " [-0.77386745]\n", + " [-0.74182428]\n", + " [-0.09146056]\n", + " [-0.83655471]\n", + " [-0.86970308]\n", + " [ 0.4025748 ]\n", + " [ 0.01829472]\n", + " [ 1.54093715]\n", + " [-0.9543623 ]\n", + " [ 0.23369796]\n", + " [-0.94776997]\n", + " [ 0.53692862]\n", + " [-0.90618744]\n", + " [ 0.0845006 ]\n", + " [ 0.30365762]\n", + " [-0.92670967]\n", + " [ 2.94938329]\n", + " [ 0.48295688]\n", + " [-0.92060875]\n", + " [ 1.36967559]\n", + " [ 0.17571516]]\n" + ] + } + ], + "source": [ + "print(gamma_activity)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualize Brain" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Plotting 64 ecog locations\n", + "Using mayavi 3d backend.\n", + "\n" + ] + } + ], + "source": [ + "subjects_dir = mne.datasets.sample.data_path() + '/subjects'\n", + "fig = plot_alignment(info, subject='sample', subjects_dir=subjects_dir,\n", + " surfaces=['pial'])\n", + "mne.viz.set_3d_view(fig, 200, 70)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "neurologic", + "language": "python", + "name": "neurologic" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index d51b7e789ed..0ba853558f7 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -14,15 +14,14 @@ # License: BSD (3-clause) import numpy as np -from sys import platform as sys_pf -# if sys_pf == 'darwin': -# import matplotlib -# matplotlib.use("Qt5Agg") import matplotlib matplotlib.use("macOSX") import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy.io import loadmat +from sklearn.preprocessing import StandardScaler + import mne from mne.viz import plot_alignment, snapshot_brain_montage @@ -47,7 +46,7 @@ print(ch_pos) # Now we make a montage stating that the ECoG contacts are in head # coordinate system (although they are in MRI). This is compensated -# by the fact that below we do not specicty a trans file so the Head<->MRI +# by the fact that below we do not specify a trans file so the Head<->MRI # transform is the identity. # montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)), # coord_frame='head') @@ -59,6 +58,40 @@ info = mne.create_info(ch_names, 1000., 'ecog').set_montage(montage) +############################################################################### +# Now that we have our electrode positions in MRI coordinates, we can load in +# our corresponding time-series data. We then compute a time-frequency +# representation of the data (i.e. 1-30, or 30-90 Hz). + +# first we'll load in the sample dataset +raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf') + +# drop bad channels +raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) + +# attach montage +raw.set_montage(montage, on_missing='warn') + +# create a 1 Epoch data structure +epoch = mne.EpochsArray(raw.get_data(start=0, stop=10000)[np.newaxis, ...], + info=raw.info) + +# perform gamma band frequency averaged over entire time period +tfr_pwr, _ = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(30, 90, 60), + n_cycles=2) +gamma_activity = tfr_pwr.data.mean(axis=(1, 2)) + +# normalize activity over all channels +# gamma_activity = StandardScaler().fit_transform(gamma_activity[:, np.newaxis]).squeeze() + +# perform low frequency activity averaged over entire time period +tfr_pwr, _ = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(8, 12, 4), + n_cycles=2) +low_activity = tfr_pwr.data.mean(axis=(1, 2)) + +# normalize activity over all channels +# low_activity = StandardScaler().fit_transform(low_activity[:, np.newaxis]).squeeze() + ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. # @@ -74,83 +107,54 @@ # Sometimes it is useful to make a scatterplot for the current figure view. # This is best accomplished with matplotlib. We can capture an image of the # current mayavi view, along with the xy position of each electrode, with the -# `snapshot_brain_montage` function. +# `snapshot_brain_montage` function. We can visualize say the gamma frequency +# of the ECoG activity on the brain using MNE functions. # We'll once again plot the surface, then take a snapshot. -fig_scatter = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, - surfaces='pial') -mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, montage) +fig_gamma = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, + surfaces='pial') +mne.viz.set_3d_view(fig_gamma, 200, 70) +xy, im = snapshot_brain_montage(fig_gamma, montage) # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) -# Define an arbitrary "activity" pattern for viz -activity = np.linspace(100, 200, xy_pts.shape[0]) - -# # This allows us to use matplotlib to create arbitrary 2d scatterplots -_, ax = plt.subplots(figsize=(10, 10)) -ax.imshow(im) -ax.scatter(*xy_pts.T, c=activity, s=200, cmap='coolwarm') -ax.set_axis_off() -# plt.show() - -############################################################################### -# Sometimes it is useful to create an animation of the ECoG activity over time. -# We can visualize say the gamma frequency of the ECoG activity on the brain -# using MNE functions. - -# first we'll load in the sample dataset -raw = mne.io.read_raw_edf(mne.datasets.misc.data_path() + '/ecog/sample_ecog.edf') - -# drop bad channels -raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) - -# attach montage -raw.set_montage(montage, on_missing='warn') - -# perform gamma band frequency -epoch = mne.EpochsArray(raw.get_data()[np.newaxis, ...], info=raw.info) -print(epoch) -# print(epoch.shape) -tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(30, 90, 60), - n_cycles=3) -print(tfr_pwr) -# Define an arbitrary "activity" pattern for viz -gamma_activity = tfr_pwr.data.mean(axis=(1, 2)) - -tfr_pwr, tfr_itc = mne.time_frequency.tfr_morlet(epoch, freqs=np.linspace(1, 30, 60), - n_cycles=3) -low_activity = tfr_pwr.data.mean(axis=(1, 2)) - - -_, ax = plt.subplots(figsize=(10, 10)) # show activity between low frequency and higher frequencies # We'll once again plot the surface, then take a snapshot. -fig_scatter = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, - surfaces='pial') -mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, montage) +fig_gamma = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, + surfaces='pial') +mne.viz.set_3d_view(fig_gamma, 200, 70) +xy, im = snapshot_brain_montage(fig_gamma, montage) +print(fig_gamma) + # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) +# show activity at higher frequencies +fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -ax.scatter(*xy_pts.T, c=gamma_activity, s=200, cmap='coolwarm') -# ax.set_title("Gamma frequency (30-90 Hz)") +sc = ax.scatter(*xy_pts.T, c=gamma_activity, s=200, cmap='viridis', + # norm=matplotlib.colors.LogNorm() + ) +ax.set_title("Gamma frequency (30-90 Hz)") +fig.colorbar(sc, ax=ax) -_, ax = plt.subplots(figsize=(10, 10)) -# show activity between low frequency and higher frequencies -# We'll once again plot the surface, then take a snapshot. -fig_scatter = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, - surfaces='pial') -mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, montage) -# Convert from a dictionary to array to plot -xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) + +# show activity between low frequency +fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -ax.scatter(*xy_pts.T, c=low_activity, s=200, cmap='coolwarm') -# ax.set_title("Low frequency (0-30 Hz)") +sc = ax.scatter(*xy_pts.T, c=low_activity, s=200, cmap='viridis', + # norm=matplotlib.colors.LogNorm() + ) +ax.set_title("Low frequency (0-30 Hz)") +fig.colorbar(sc, ax=ax) + +# create an axes on the right side of ax. The width of cax will be 5% +# of ax and the padding between cax and ax will be fixed at 0.05 inch. +# divider = make_axes_locatable(ax) +# cax = divider.append_axes("right", size="5%", pad=0.05) +# plt.colorbar(sc, ax=cax) plt.show() From 1a9b7c8ea63dc00163728f484c781abecf5e3b00 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 12:10:05 -0400 Subject: [PATCH 07/27] Fix on missing. --- tutorials/misc/plot_ecog.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index b158f3d98bf..709d33966be 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -98,7 +98,7 @@ def _from_tsv(fname): raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) # attach montage -raw.set_montage(montage, on_missing='warn') +raw.set_montage(montage, on_missing='ignore') # create a 1 Epoch data structure epoch = mne.EpochsArray(raw.get_data()[np.newaxis, ...], From 02d9a517c6fc42c288f0fd4f902d889d29c9e1c9 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 15:37:43 -0400 Subject: [PATCH 08/27] Updating misc path. --- mne/datasets/utils.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index d27deb21332..7731d560963 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -238,7 +238,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, path = _get_path(path, key, name) # To update the testing or misc dataset, push commits, then make a new # release on GitHub. Then update the "releases" variable: - releases = dict(testing='0.85', misc='0.5') + releases = dict(testing='0.85', misc='0.6') # And also update the "md5_hashes['testing']" variable below. # To update any other dataset, update the data archive itself (upload @@ -254,9 +254,8 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, bst_resting='https://osf.io/m7bd3/download?version=3'), fake='https://github.com/mne-tools/mne-testing-data/raw/master/' 'datasets/foo.tgz', - # misc='https://codeload.github.com/mne-tools/mne-misc-data/' - # 'tar.gz/%s' % releases['misc'], - misc='https://github.com/adam2392/mne-misc-data/tarball/seeg', + misc='https://codeload.github.com/mne-tools/mne-misc-data/' + 'tar.gz/%s' % releases['misc'], sample='https://osf.io/86qa2/download?version=5', somato='https://osf.io/tp4sg/download?version=6', spm='https://osf.io/je4s8/download?version=2', @@ -278,8 +277,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, archive_names = dict( fieldtrip_cmc='SubjectCMC.zip', kiloword='MNE-kiloword-data.tar.gz', - # misc='mne-misc-data-%s.tar.gz' % releases['misc'], - misc='mne-misc-data.tar.gz', + misc='mne-misc-data-%s.tar.gz' % releases['misc'], mtrf='mTRF_1.5.zip', multimodal='MNE-multimodal-data.tar.gz', fnirs_motor='MNE-fNIRS-motor-data.tgz', @@ -295,8 +293,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, # original folder names that get extracted (only needed if the # archive does not extract the right folder name; e.g., usually GitHub) folder_origs = dict( # not listed means None (no need to move) - # misc='mne-misc-data-%s' % releases['misc'], - misc='adam2392-mne-misc-data-75968e3', + misc='mne-misc-data-%s' % releases['misc'], testing='mne-testing-data-%s' % releases['testing'], ) # finally, where we want them to extract to (only needed if the folder name @@ -321,7 +318,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, bst_raw='fa2efaaec3f3d462b319bc24898f440c', bst_resting='70fc7bf9c3b97c4f2eab6260ee4a0430'), fake='3194e9f7b46039bb050a74f3e1ae9908', - misc='a82a6888d4266d3c6d5aacc8f305b5b8', + misc='e00808c3b05123059e2cf49ff276e919', sample='12b75d1cb7df9dfb4ad73ed82f61094f', somato='ea825966c0a1e9b2f84e3826c5500161', spm='9f43f67150e3b694b523a21eb929ea75', From f03b2b7e9e8317914571aefb14e1afd1344b2c7b Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 17:13:33 -0400 Subject: [PATCH 09/27] Updating change log. --- doc/changes/latest.inc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 214e1648e6e..63cd950bac6 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -61,6 +61,8 @@ Changelog - Add ``plot`` option to :meth:`mne.viz.plot_filter` allowing selection of which filter properties are plotted and added option for user to supply ``axes`` by `Robert Luke`_ +- Add ECoG misc EDF datast to the :ref:`tut-plot_ecog` to show snapshots of time-frequency acitivity by `Adam Li`_ + Bug ~~~ From ef656df7f871b4ff4da9dd1391ae97a94c7be5e3 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Tue, 12 May 2020 20:58:58 -0400 Subject: [PATCH 10/27] Adding name to plot ecog example. --- doc/changes/latest.inc | 2 +- tutorials/misc/plot_ecog.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 7e6bba55ff9..0b3399ac364 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -67,7 +67,7 @@ Changelog - Add ``sources`` and ``detectors`` options for fNIRS use of :meth:`mne.viz.plot_alignment` allowing plotting of optode locations in addition to channel midpoint ``channels`` and ``path`` between fNIRS optodes by `Kyle Mathewson`_ -- Add ECoG misc EDF datast to the :ref:`tut-plot_ecog` to show snapshots of time-frequency acitivity by `Adam Li`_ +- Add ECoG misc EDF dataset to the :ref:`tut_working_with_ecog` tutorial to show snapshots of time-frequency activity by `Adam Li`_ Bug ~~~ diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 709d33966be..62adf8d84fe 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -1,4 +1,6 @@ """ +.. _tut_working_with_ecog: + ====================== Working with ECoG data ====================== From 3b7c53ad7f6a7821a249d9c9b1dfe455afbc791a Mon Sep 17 00:00:00 2001 From: Adam Li Date: Wed, 13 May 2020 12:13:14 -0400 Subject: [PATCH 11/27] Adjust agram comments. --- tutorials/misc/plot_ecog.py | 41 +++++++------------------------------ 1 file changed, 7 insertions(+), 34 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 62adf8d84fe..8dd2f4f7cf1 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -11,11 +11,12 @@ """ # Authors: Eric Larson # Chris Holdgraf -# Edited: Adam Li +# Adam Li # # License: BSD (3-clause) -from collections import OrderedDict +import collections +import pandas as pd import numpy as np import matplotlib.pyplot as plt @@ -34,40 +35,12 @@ # a :class:`mne.channels.DigMontage` class. # First, define a helper function to read in .tsv file. - -def _from_tsv(fname): - """Read a tsv file into an OrderedDict. - - - Parameters - ---------- - fname : str - Path to the file being loaded. - - Returns - ------- - data_dict : collections.OrderedDict - Keys are the column names, and values are the column data. - - """ - data = np.loadtxt(fname, dtype=str, delimiter='\t', - comments=None, encoding='utf-8') - column_names = data[0, :] - info = data[1:, :] - data_dict = OrderedDict() - dtypes = [str] * info.shape[1] - for i, name in enumerate(column_names): - data_dict[name] = info[:, i].astype(dtypes[i]).tolist() - return data_dict - - # read in the electrode coordinates file # in this tutorial, these are assumed to be in meters -elec_tsv = _from_tsv(misc_path + '/ecog/sample_ecog_electrodes.tsv') -ch_names = elec_tsv['name'] -ch_coords = np.vstack((elec_tsv['x'], - elec_tsv['y'], - elec_tsv['z'])).T.astype(float) +elec_df = pd.read_csv(misc_path + '/ecog/sample_ecog_electrodes.tsv', + sep='\t', header=0, index_col=None) +ch_names = elec_df['name'].tolist() +ch_coords = elec_df[['x', 'y', 'z']].to_numpy(dtype=float) ch_pos = dict(zip(ch_names, ch_coords)) # create montage from channel coordinates in the 'head' coordinate frame From 4bf98697d04ec5e52c3599e551a265b32ef72f72 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Wed, 13 May 2020 12:13:24 -0400 Subject: [PATCH 12/27] Adjust agram comments. --- tutorials/misc/plot_ecog.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 8dd2f4f7cf1..e5e550cf809 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -15,7 +15,6 @@ # # License: BSD (3-clause) -import collections import pandas as pd import numpy as np import matplotlib.pyplot as plt From cffef022bcde62eb4c3b25eed50ae0cfb00a4994 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Wed, 13 May 2020 16:39:36 -0400 Subject: [PATCH 13/27] Cleanup. --- tutorials/misc/plot_ecog.py | 49 +++++++++++++++---------------------- 1 file changed, 20 insertions(+), 29 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index e5e550cf809..5aa9e71006f 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -20,7 +20,6 @@ import matplotlib.pyplot as plt import mne -from mne.time_frequency import tfr_morlet from mne.viz import plot_alignment, snapshot_brain_montage print(__doc__) @@ -42,16 +41,13 @@ ch_coords = elec_df[['x', 'y', 'z']].to_numpy(dtype=float) ch_pos = dict(zip(ch_names, ch_coords)) -# create montage from channel coordinates in the 'head' coordinate frame -montage = mne.channels.make_dig_montage(ch_pos, - coord_frame='head') - # Now we make a montage stating that the ECoG contacts are in head # coordinate system (although they are in MRI). This is compensated # by the fact that below we do not specify a trans file so the Head<->MRI # transform is the identity. -# montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)), -# coord_frame='head') + +# create montage from channel coordinates in the 'head' coordinate frame +montage = mne.channels.make_dig_montage(ch_pos, coord_frame='head') print('Created %s channel positions' % len(ch_names)) ############################################################################### @@ -63,30 +59,22 @@ ############################################################################### # Now that we have our electrode positions in MRI coordinates, we can load in # our corresponding time-series data. We then compute a time-frequency -# representation of the data (i.e. 1-30, or 30-90 Hz). +# representation of the data (i.e. 30-90 Hz). # first we'll load in the sample dataset raw = mne.io.read_raw_edf(misc_path + '/ecog/sample_ecog.edf') # drop bad channels raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) +raw.load_data() +raw.drop_channels(raw.info['bads']) # attach montage -raw.set_montage(montage, on_missing='ignore') - -# create a 1 Epoch data structure -epoch = mne.EpochsArray(raw.get_data()[np.newaxis, ...], - info=raw.info) - -# perform gamma band frequency averaged over entire time period -tfr_pwr, _ = tfr_morlet(epoch, freqs=np.linspace(30, 90, 60), - n_cycles=2) -gamma_activity = tfr_pwr.data.mean(axis=(1, 2)) +raw.set_montage(montage) -# perform low frequency activity averaged over entire time period -tfr_pwr, _ = tfr_morlet(epoch, freqs=np.linspace(8, 12, 4), - n_cycles=2) -low_activity = tfr_pwr.data.mean(axis=(1, 2)) +# compute gamma and alpha band activity +gamma_activity = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) +alpha_activity = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. @@ -115,26 +103,29 @@ # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) -vmin, vmax = np.percentile(gamma_activity, [10, 90]) +# colormap to view spectral power +cmap = 'viridis' # show activity at higher frequencies +vmin, vmax = np.percentile(gamma_activity, [10, 90]) + fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() sc = ax.scatter(*xy_pts.T, c=gamma_activity, s=200, - cmap='viridis', vmin=vmin, vmax=vmax) + cmap=cmap, vmin=vmin, vmax=vmax) ax.set_title("Gamma frequency (30-90 Hz)") fig.colorbar(sc, ax=ax) -vmin, vmax = np.percentile(low_activity, [10, 90]) - # show activity between low frequency +vmin, vmax = np.percentile(alpha_activity, [10, 90]) + fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -sc = ax.scatter(*xy_pts.T, c=low_activity, s=200, - cmap='viridis', vmin=vmin, vmax=vmax) -ax.set_title("Low frequency (0-30 Hz)") +sc = ax.scatter(*xy_pts.T, c=alpha_activity, s=200, + cmap=cmap, vmin=vmin, vmax=vmax) +ax.set_title("Alpha frequency (8-12 Hz)") fig.colorbar(sc, ax=ax) plt.show() From be1aa906e3b2a8e3a90fd411879d7d3242f6b554 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Wed, 13 May 2020 16:40:24 -0400 Subject: [PATCH 14/27] Cleanup. --- tutorials/misc/plot_ecog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 5aa9e71006f..e1acfaf32c7 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -31,10 +31,10 @@ ############################################################################### # Let's load some ECoG electrode locations and names, and turn them into # a :class:`mne.channels.DigMontage` class. -# First, define a helper function to read in .tsv file. +# First, use pandas to read in the .tsv file. # read in the electrode coordinates file -# in this tutorial, these are assumed to be in meters +# in this tutorial, coordinates are assumed to be in meters elec_df = pd.read_csv(misc_path + '/ecog/sample_ecog_electrodes.tsv', sep='\t', header=0, index_col=None) ch_names = elec_df['name'].tolist() From afcc5b06dbdf00aa641a5d4eb4fb428bff783a81 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Thu, 14 May 2020 12:02:38 +0200 Subject: [PATCH 15/27] pass on tuto --- tutorials/misc/plot_ecog.py | 51 ++++++++++++++++--------------------- 1 file changed, 22 insertions(+), 29 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index e1acfaf32c7..efd533e6c51 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -33,33 +33,24 @@ # a :class:`mne.channels.DigMontage` class. # First, use pandas to read in the .tsv file. -# read in the electrode coordinates file -# in this tutorial, coordinates are assumed to be in meters +# In this tutorial, the electrode coordinates are assumed to be in meters elec_df = pd.read_csv(misc_path + '/ecog/sample_ecog_electrodes.tsv', sep='\t', header=0, index_col=None) ch_names = elec_df['name'].tolist() ch_coords = elec_df[['x', 'y', 'z']].to_numpy(dtype=float) ch_pos = dict(zip(ch_names, ch_coords)) -# Now we make a montage stating that the ECoG contacts are in head -# coordinate system (although they are in MRI). This is compensated -# by the fact that below we do not specify a trans file so the Head<->MRI -# transform is the identity. +# Now we make a :class:`mne.channels.DigMontage` stating that the ECoG +# contacts are in head coordinate system (although they are in MRI). This is +# compensated below by the fact that we do not specify a trans file so the +# Head<->MRI transform is the identity. -# create montage from channel coordinates in the 'head' coordinate frame montage = mne.channels.make_dig_montage(ch_pos, coord_frame='head') print('Created %s channel positions' % len(ch_names)) ############################################################################### -# Now that we have our electrode positions in MRI coordinates, we can create -# our measurement info structure. - -info = mne.create_info(ch_names, 1000., 'ecog').set_montage(montage) - -############################################################################### -# Now that we have our electrode positions in MRI coordinates, we can load in -# our corresponding time-series data. We then compute a time-frequency -# representation of the data (i.e. 30-90 Hz). +# Now that we have our montage, we can load in our corresponding +# time-series data and set the montage to the raw data. # first we'll load in the sample dataset raw = mne.io.read_raw_edf(misc_path + '/ecog/sample_ecog.edf') @@ -72,9 +63,11 @@ # attach montage raw.set_montage(montage) -# compute gamma and alpha band activity -gamma_activity = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) -alpha_activity = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) +############################################################################### +# We then compute the signal power in certain frequency bands (e.g. 30-90 Hz). +# We compute the power in gamma and alpha bands. +gamma_power = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) +alpha_power = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. @@ -83,7 +76,7 @@ # do not align to the cortical surface perfectly. subjects_dir = sample_path + '/subjects' -fig = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, +fig = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, surfaces=['pial']) mne.viz.set_3d_view(fig, 200, 70) @@ -91,8 +84,8 @@ # Sometimes it is useful to make a scatterplot for the current figure view. # This is best accomplished with matplotlib. We can capture an image of the # current mayavi view, along with the xy position of each electrode, with the -# `snapshot_brain_montage` function. We can visualize say the gamma frequency -# of the ECoG activity on the brain using MNE functions. +# `snapshot_brain_montage` function. We can then visualize for example the +# gamma power on the brain using MNE and matplotlib functions. # We'll once again plot the surface, then take a snapshot. fig_scatter = plot_alignment(raw.info, subject='sample', @@ -101,29 +94,29 @@ xy, im = snapshot_brain_montage(fig_scatter, montage) # Convert from a dictionary to array to plot -xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) +xy_pts = np.vstack([xy[ch] for ch in raw.info['ch_names']]) # colormap to view spectral power cmap = 'viridis' -# show activity at higher frequencies -vmin, vmax = np.percentile(gamma_activity, [10, 90]) +# show power at higher frequencies +vmin, vmax = np.percentile(gamma_power, [10, 90]) fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -sc = ax.scatter(*xy_pts.T, c=gamma_activity, s=200, +sc = ax.scatter(*xy_pts.T, c=gamma_power, s=200, cmap=cmap, vmin=vmin, vmax=vmax) ax.set_title("Gamma frequency (30-90 Hz)") fig.colorbar(sc, ax=ax) -# show activity between low frequency -vmin, vmax = np.percentile(alpha_activity, [10, 90]) +# show power between low frequency +vmin, vmax = np.percentile(alpha_power, [10, 90]) fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -sc = ax.scatter(*xy_pts.T, c=alpha_activity, s=200, +sc = ax.scatter(*xy_pts.T, c=alpha_power, s=200, cmap=cmap, vmin=vmin, vmax=vmax) ax.set_title("Alpha frequency (8-12 Hz)") fig.colorbar(sc, ax=ax) From 140b73291aa72f1b174bb208ccd0deea02dce013 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 14 May 2020 14:55:36 -0400 Subject: [PATCH 16/27] Adding animation back in. --- tutorials/misc/plot_ecog.py | 101 ++++++++++++++++++++++++------------ 1 file changed, 68 insertions(+), 33 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index efd533e6c51..611ff3cfb8d 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -18,6 +18,7 @@ import pandas as pd import numpy as np import matplotlib.pyplot as plt +import matplotlib.animation as animation import mne from mne.viz import plot_alignment, snapshot_brain_montage @@ -40,6 +41,7 @@ ch_coords = elec_df[['x', 'y', 'z']].to_numpy(dtype=float) ch_pos = dict(zip(ch_names, ch_coords)) +############################################################################### # Now we make a :class:`mne.channels.DigMontage` stating that the ECoG # contacts are in head coordinate system (although they are in MRI). This is # compensated below by the fact that we do not specify a trans file so the @@ -63,35 +65,36 @@ # attach montage raw.set_montage(montage) -############################################################################### -# We then compute the signal power in certain frequency bands (e.g. 30-90 Hz). -# We compute the power in gamma and alpha bands. -gamma_power = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) -alpha_power = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) - ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. # # .. note:: These are not real electrodes for this subject, so they # do not align to the cortical surface perfectly. +# +# Sometimes it is useful to make a scatterplot for the current figure view. +# This is best accomplished with matplotlib. We can capture an image of the +# current mayavi view, along with the xy position of each electrode, with the +# `snapshot_brain_montage` function. We can then visualize for example the +# gamma power on the brain using MNE and matplotlib functions. +# Here, we'll use :func:`~mne.viz.snapshot_brain_montage` to save the plot +# as image data, so that later in the tutorial we can plot data on top of it. subjects_dir = sample_path + '/subjects' fig = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, surfaces=['pial']) mne.viz.set_3d_view(fig, 200, 70) +xy, im = snapshot_brain_montage(fig, montage) + ############################################################################### -# Sometimes it is useful to make a scatterplot for the current figure view. -# This is best accomplished with matplotlib. We can capture an image of the -# current mayavi view, along with the xy position of each electrode, with the -# `snapshot_brain_montage` function. We can then visualize for example the -# gamma power on the brain using MNE and matplotlib functions. +# We then compute the signal power in certain frequency bands (e.g. 30-90 Hz). +# We compute the power in gamma and alpha bands. +gamma_power = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) +alpha_power = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) -# We'll once again plot the surface, then take a snapshot. -fig_scatter = plot_alignment(raw.info, subject='sample', - subjects_dir=subjects_dir, surfaces='pial') -mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, montage) +############################################################################### +# Now let's use matplotlib to overplot frequency band power onto the electrodes +# which can be plotted on top of the brain from `snapshot_brain_montage`. # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in raw.info['ch_names']]) @@ -99,26 +102,58 @@ # colormap to view spectral power cmap = 'viridis' -# show power at higher frequencies -vmin, vmax = np.percentile(gamma_power, [10, 90]) - -fig, ax = plt.subplots(figsize=(10, 10)) -ax.imshow(im) -ax.set_axis_off() -sc = ax.scatter(*xy_pts.T, c=gamma_power, s=200, - cmap=cmap, vmin=vmin, vmax=vmax) -ax.set_title("Gamma frequency (30-90 Hz)") +# create a 1x2 subplot showing a snapshot of the average power +# in gamma and alpha frequency bands +fig, axs = plt.subplots(1, 2, figsize=(20, 10)) +_gamma_alpha_power = np.concatenate((gamma_power, alpha_power)).flatten() +vmin, vmax = np.percentile(_gamma_alpha_power, [10, 90]) +for ax, band_power, band in zip(axs, + [gamma_power, alpha_power], + ['Gamma', 'Alpha']): + ax.imshow(im) + ax.set_axis_off() + sc = ax.scatter(*xy_pts.T, c=band_power, s=200, + cmap=cmap, vmin=vmin, vmax=vmax) + ax.set_title(f'{band} frequency', size='large') fig.colorbar(sc, ax=ax) -# show power between low frequency -vmin, vmax = np.percentile(alpha_power, [10, 90]) +############################################################################### +# Say we want to visualize the evolution of the power in the gamma band, +# instead of just plotting the average. We can use +# `matplotlib.animation.FuncAnimation` to create an animation and apply this +# to the brain figure. + +# create an initialization and animation function +# to pass to FuncAnimation +def init(): + """Create an empty frame.""" + return paths, + + +def animate(i, activity): + """Animate the plot.""" + paths.set_array(activity[:, i]) + return paths, + + +# create an Epoch to use morlet wavelet transform +events = [[raw.first_samp, 0, 1]] +epochs = mne.Epochs(raw, events, tmin=0, tmax=raw.times[-1], + preload=True) + +# compute time series of the gamma frequency band +tfr_pwr, _ = mne.time_frequency.tfr_morlet( + epochs, freqs=np.linspace(30, 90), n_cycles=7) +gamma_activity = tfr_pwr.data.mean(axis=1) + +# create the figure and apply the animation of the +# gamma frequency band activity fig, ax = plt.subplots(figsize=(10, 10)) ax.imshow(im) ax.set_axis_off() -sc = ax.scatter(*xy_pts.T, c=alpha_power, s=200, - cmap=cmap, vmin=vmin, vmax=vmax) -ax.set_title("Alpha frequency (8-12 Hz)") -fig.colorbar(sc, ax=ax) - -plt.show() +paths = ax.scatter(*xy_pts.T, c=np.zeros(len(xy_pts)), s=200, + cmap=cmap, vmin=vmin, vmax=vmax) +anim = animation.FuncAnimation(fig, animate, init_func=init, + fargs=(gamma_activity,), frames=500, + interval=20, blit=True) From 3d87f1cd4e32fbe1c7d68305ba4254be4439ede7 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 14 May 2020 16:06:29 -0400 Subject: [PATCH 17/27] Fixing error in epoch creation. --- tutorials/misc/plot_ecog.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 611ff3cfb8d..9fc077a57bf 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -140,7 +140,7 @@ def animate(i, activity): # create an Epoch to use morlet wavelet transform events = [[raw.first_samp, 0, 1]] epochs = mne.Epochs(raw, events, tmin=0, tmax=raw.times[-1], - preload=True) + baseline=None, preload=True) # compute time series of the gamma frequency band tfr_pwr, _ = mne.time_frequency.tfr_morlet( From 33a533a6e1ec64c9d2cb5727262d133a0dd820f4 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Thu, 14 May 2020 16:56:31 -0400 Subject: [PATCH 18/27] Longer frames. --- tutorials/misc/plot_ecog.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 9fc077a57bf..786e5b94a5f 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -154,6 +154,10 @@ def animate(i, activity): ax.set_axis_off() paths = ax.scatter(*xy_pts.T, c=np.zeros(len(xy_pts)), s=200, cmap=cmap, vmin=vmin, vmax=vmax) +fig.colorbar(paths, ax=ax) +ax.set_title('Gamma frequency over time (Morlet wavelet)', + size='large') anim = animation.FuncAnimation(fig, animate, init_func=init, - fargs=(gamma_activity,), frames=500, + fargs=(gamma_activity,), + frames=gamma_activity.shape[1], interval=20, blit=True) From a494838c60a9a7e189c764d357d7a4c0e0ef975c Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 14 May 2020 21:14:16 -0400 Subject: [PATCH 19/27] FIX: Only one power calc --- tutorials/misc/plot_ecog.py | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 786e5b94a5f..73a48f24dd2 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -61,6 +61,7 @@ raw.info['bads'].extend([ch for ch in raw.ch_names if ch not in ch_names]) raw.load_data() raw.drop_channels(raw.info['bads']) +raw.crop(0, 2) # just process 2 sec of data for speed # attach montage raw.set_montage(montage) @@ -89,9 +90,13 @@ ############################################################################### # We then compute the signal power in certain frequency bands (e.g. 30-90 Hz). # We compute the power in gamma and alpha bands. -gamma_power = np.sum(raw.copy().filter(30, 90).get_data() ** 2, axis=1) -alpha_power = np.sum(raw.copy().filter(8, 12).get_data() ** 2, axis=1) - +gamma_power_t = raw.copy().filter(30, 90).apply_hilbert( + envelope=True).get_data() +alpha_power_t = raw.copy().filter(8, 12).apply_hilbert( + envelope=True).get_data() +gamma_power = gamma_power_t.mean(axis=-1) +alpha_power = alpha_power_t.mean(axis=-1) +# these have edge artifcat ############################################################################### # Now let's use matplotlib to overplot frequency band power onto the electrodes # which can be plotted on top of the brain from `snapshot_brain_montage`. @@ -137,16 +142,6 @@ def animate(i, activity): return paths, -# create an Epoch to use morlet wavelet transform -events = [[raw.first_samp, 0, 1]] -epochs = mne.Epochs(raw, events, tmin=0, tmax=raw.times[-1], - baseline=None, preload=True) - -# compute time series of the gamma frequency band -tfr_pwr, _ = mne.time_frequency.tfr_morlet( - epochs, freqs=np.linspace(30, 90), n_cycles=7) -gamma_activity = tfr_pwr.data.mean(axis=1) - # create the figure and apply the animation of the # gamma frequency band activity fig, ax = plt.subplots(figsize=(10, 10)) @@ -157,7 +152,10 @@ def animate(i, activity): fig.colorbar(paths, ax=ax) ax.set_title('Gamma frequency over time (Morlet wavelet)', size='large') + +# avoid edge artifacts and decimate, showing just a short chunk +show_power = gamma_power_t[:, 100:-1700:2] anim = animation.FuncAnimation(fig, animate, init_func=init, - fargs=(gamma_activity,), - frames=gamma_activity.shape[1], - interval=20, blit=True) + fargs=(show_power,), + frames=show_power.shape[1], + interval=0.2, blit=True) From 70df9eaa180a82f298cb12f23f0a38f982c78258 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 14 May 2020 21:29:22 -0400 Subject: [PATCH 20/27] FIX: Interval --- tutorials/misc/plot_ecog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 73a48f24dd2..666a615dd97 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -96,7 +96,7 @@ envelope=True).get_data() gamma_power = gamma_power_t.mean(axis=-1) alpha_power = alpha_power_t.mean(axis=-1) -# these have edge artifcat + ############################################################################### # Now let's use matplotlib to overplot frequency band power onto the electrodes # which can be plotted on top of the brain from `snapshot_brain_montage`. @@ -158,4 +158,4 @@ def animate(i, activity): anim = animation.FuncAnimation(fig, animate, init_func=init, fargs=(show_power,), frames=show_power.shape[1], - interval=0.2, blit=True) + interval=100, blit=True) From f1d6530fe335de332d12d79ae8c95a382b2d9075 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 14 May 2020 21:33:41 -0400 Subject: [PATCH 21/27] FIX: Thumb --- tutorials/misc/plot_ecog.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 666a615dd97..20ff0130d3c 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -150,9 +150,10 @@ def animate(i, activity): paths = ax.scatter(*xy_pts.T, c=np.zeros(len(xy_pts)), s=200, cmap=cmap, vmin=vmin, vmax=vmax) fig.colorbar(paths, ax=ax) -ax.set_title('Gamma frequency over time (Morlet wavelet)', +ax.set_title('Gamma frequency over time (Hilbert transform)', size='large') +# sphinx_gallery_thumbnail_number = 3 # avoid edge artifacts and decimate, showing just a short chunk show_power = gamma_power_t[:, 100:-1700:2] anim = animation.FuncAnimation(fig, animate, init_func=init, From 9a050a01647b50c9e899bcf13d9adbc1d678c3c4 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:17 -0400 Subject: [PATCH 22/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 20ff0130d3c..2948f278b72 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -68,17 +68,12 @@ ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. +# We'll use :func:`~mne.viz.snapshot_brain_montage` to save the plot as image +# data (along with xy positions of each electrode in the image), so that later +# we can plot frequency band power on top of it. # # .. note:: These are not real electrodes for this subject, so they # do not align to the cortical surface perfectly. -# -# Sometimes it is useful to make a scatterplot for the current figure view. -# This is best accomplished with matplotlib. We can capture an image of the -# current mayavi view, along with the xy position of each electrode, with the -# `snapshot_brain_montage` function. We can then visualize for example the -# gamma power on the brain using MNE and matplotlib functions. -# Here, we'll use :func:`~mne.viz.snapshot_brain_montage` to save the plot -# as image data, so that later in the tutorial we can plot data on top of it. subjects_dir = sample_path + '/subjects' fig = plot_alignment(raw.info, subject='sample', subjects_dir=subjects_dir, From b71595c63e65df474608a46f744880e1213bb049 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:26 -0400 Subject: [PATCH 23/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 2948f278b72..121fd70d809 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -94,7 +94,8 @@ ############################################################################### # Now let's use matplotlib to overplot frequency band power onto the electrodes -# which can be plotted on top of the brain from `snapshot_brain_montage`. +# which can be plotted on top of the brain from +# :func:`~mne.viz.snapshot_brain_montage`. # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in raw.info['ch_names']]) From 9b154bdfca46521e4372e1ca2ed0e36e4a7aed1f Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:36 -0400 Subject: [PATCH 24/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 121fd70d809..d1995b5a9fb 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -83,8 +83,8 @@ xy, im = snapshot_brain_montage(fig, montage) ############################################################################### -# We then compute the signal power in certain frequency bands (e.g. 30-90 Hz). -# We compute the power in gamma and alpha bands. +# Next, we'll compute the signal power in the gamma (30-90 Hz) and alpha +# (8-12 Hz) bands. gamma_power_t = raw.copy().filter(30, 90).apply_hilbert( envelope=True).get_data() alpha_power_t = raw.copy().filter(8, 12).apply_hilbert( From ea7441c53f841c4dff3667143c583030250511c7 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:44 -0400 Subject: [PATCH 25/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index d1995b5a9fb..f5b6c3cd54c 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -106,6 +106,7 @@ # create a 1x2 subplot showing a snapshot of the average power # in gamma and alpha frequency bands fig, axs = plt.subplots(1, 2, figsize=(20, 10)) +# choose a colormap range wide enough for both frequency bands _gamma_alpha_power = np.concatenate((gamma_power, alpha_power)).flatten() vmin, vmax = np.percentile(_gamma_alpha_power, [10, 90]) for ax, band_power, band in zip(axs, From 1cbca05512fd0ef0bf59bb2fec528be5f37524e3 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:51 -0400 Subject: [PATCH 26/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index f5b6c3cd54c..b7269e052b3 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -116,8 +116,8 @@ ax.set_axis_off() sc = ax.scatter(*xy_pts.T, c=band_power, s=200, cmap=cmap, vmin=vmin, vmax=vmax) - ax.set_title(f'{band} frequency', size='large') -fig.colorbar(sc, ax=ax) + ax.set_title(f'{band} band power', size='x-large') +fig.colorbar(sc, ax=axs) ############################################################################### # Say we want to visualize the evolution of the power in the gamma band, From 6e147c6fa4105951ea87e11181047b898f1321d9 Mon Sep 17 00:00:00 2001 From: Adam Li Date: Fri, 15 May 2020 13:05:58 -0400 Subject: [PATCH 27/27] Update tutorials/misc/plot_ecog.py Co-authored-by: Daniel McCloy --- tutorials/misc/plot_ecog.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index b7269e052b3..9f85f26cc7d 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -103,8 +103,7 @@ # colormap to view spectral power cmap = 'viridis' -# create a 1x2 subplot showing a snapshot of the average power -# in gamma and alpha frequency bands +# Create a 1x2 figure showing the average power in gamma and alpha bands. fig, axs = plt.subplots(1, 2, figsize=(20, 10)) # choose a colormap range wide enough for both frequency bands _gamma_alpha_power = np.concatenate((gamma_power, alpha_power)).flatten()