From 244e2ccde57cd5c2ac3de3c863a4bc71f730abad Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Tue, 1 Oct 2019 14:43:31 -0400 Subject: [PATCH 01/24] source-snr-edit --- mne/source_estimate.py | 80 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 79 insertions(+), 1 deletion(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 73eefd2e769..9b4d6d08da4 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -25,7 +25,8 @@ plot_volume_source_estimates) from .io.base import ToDataFrameMixin, TimeMixin from .externals.h5io import read_hdf5, write_hdf5 - +#from mne.forward import convert_forward_solution +#from mne.minimum_norm.inverse import _prepare_forward def _read_stc(filename): """Aux Function.""" @@ -2744,3 +2745,80 @@ def extract_label_time_course(stcs, labels, src, mode='mean_flip', label_tc = label_tc[0] return label_tc + + +@verbose +def estimate_snr(self, info, fwd, cov, verbose=None): + """Compute time-varying SNR + + This function should only be used with source estimates with units + nA (i.e., MNE-like solutions, *not* dSPM or sLORETA). + + Reference: + Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., + Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + Mapping the Signal-To-Noise-Ratios of Cortical Sources in + Magnetoencephalography and Electroencephalography. + Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 + + Parameters + ---------- + info : instance of io.meas_info.Info + The measurement info. + fwd : instance of Forward + The forward solution used to create the source estimate. + cov : instance of Covariance + The noise covariance used to estimate the resting cortical + activations. Should be an evoked covariance, not empty room. + + Returns + ------- + snr_stc : instance of SourceEstimate + The source estimate with the SNR computed. + + Notes + ----- + We define the SNR in decibels for each source location at each + time point as: + + .. math:: + + {\\rm SNR} = 10\\log_10[\\frac{a^2}{N}\\sum_k\\frac{b_k^2}{s_k^2}] + + where :math:`\\b_k` is the signal on sensor :math:`k` provided by the + forward model for a source with unit amplitude, :math:`a` is the + source amplitude, :math:`N` is the number of sensors, and + :math:`s_k^2` is the noise variance on sensor :math:`k`. + """ + from .forward import convert_forward_solution + from .minimum_norm.inverse import _prepare_forward + + fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) + print("done convert forward solution\n") + + + + #def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, + # use_cps, exp, limit_depth_chs, combine_xyz, + # allow_fixed_depth, limit): + +# return (forward, info_picked, gain, depth_prior, orient_prior, source_std, +# trace_GRGT, noise_cov, whitener) + + _, _, _, _, _, A, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, + use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', + allow_fixed_depth=True, limit=None) + + print("done _prepare_forward\n") + N = cov['dim'] + print ("N = %d sensors\n" % N) + b_k2 = (A * A).T + s_k2 = np.diag(cov['data']) + + print (np.shape(b_k2)) + print (np.shape(s_k2)) + scaling = (1. / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] +# print (np.shape(scaling)) + + print("End of SNR Func\n") + \ No newline at end of file From 908e7b9d9313a1c41e45d006bafcf80be3d05072 Mon Sep 17 00:00:00 2001 From: Mainak Jas Date: Tue, 1 Oct 2019 15:02:49 -0400 Subject: [PATCH 02/24] FIX move estimate_snr to be a method --- mne/source_estimate.py | 152 +++++++++--------- .../source-modeling/plot_mne_solutions.py | 4 + 2 files changed, 79 insertions(+), 77 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 9b4d6d08da4..bc6927f4149 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1099,6 +1099,81 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False): return stcs + @verbose + def estimate_snr(self, info, fwd, cov, verbose=None): + """Compute time-varying SNR + + This function should only be used with source estimates with units + nA (i.e., MNE-like solutions, *not* dSPM or sLORETA). + + Reference: + Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., + Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + Mapping the Signal-To-Noise-Ratios of Cortical Sources in + Magnetoencephalography and Electroencephalography. + Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 + + Parameters + ---------- + info : instance of io.meas_info.Info + The measurement info. + fwd : instance of Forward + The forward solution used to create the source estimate. + cov : instance of Covariance + The noise covariance used to estimate the resting cortical + activations. Should be an evoked covariance, not empty room. + + Returns + ------- + snr_stc : instance of SourceEstimate + The source estimate with the SNR computed. + + Notes + ----- + We define the SNR in decibels for each source location at each + time point as: + + .. math:: + + {\\rm SNR} = 10\\log_10[\\frac{a^2}{N}\\sum_k\\frac{b_k^2}{s_k^2}] + + where :math:`\\b_k` is the signal on sensor :math:`k` provided by the + forward model for a source with unit amplitude, :math:`a` is the + source amplitude, :math:`N` is the number of sensors, and + :math:`s_k^2` is the noise variance on sensor :math:`k`. + """ + from .forward import convert_forward_solution + from .minimum_norm.inverse import _prepare_forward + + fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) + print("done convert forward solution\n") + + + + #def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, + # use_cps, exp, limit_depth_chs, combine_xyz, + # allow_fixed_depth, limit): + + # return (forward, info_picked, gain, depth_prior, orient_prior, source_std, + # trace_GRGT, noise_cov, whitener) + + _, _, _, _, _, A, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, + use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', + allow_fixed_depth=True, limit=None) + + print("done _prepare_forward\n") + N = cov['dim'] + print ("N = %d sensors\n" % N) + b_k2 = (A * A).T + s_k2 = np.diag(cov['data']) + + print (np.shape(b_k2)) + print (np.shape(s_k2)) + scaling = (1. / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] + # print (np.shape(scaling)) + + print("End of SNR Func\n") + def _center_of_mass(vertices, values, hemi, surf, subject, subjects_dir, restrict_vertices): @@ -2745,80 +2820,3 @@ def extract_label_time_course(stcs, labels, src, mode='mean_flip', label_tc = label_tc[0] return label_tc - - -@verbose -def estimate_snr(self, info, fwd, cov, verbose=None): - """Compute time-varying SNR - - This function should only be used with source estimates with units - nA (i.e., MNE-like solutions, *not* dSPM or sLORETA). - - Reference: - Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., - Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). - Mapping the Signal-To-Noise-Ratios of Cortical Sources in - Magnetoencephalography and Electroencephalography. - Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 - - Parameters - ---------- - info : instance of io.meas_info.Info - The measurement info. - fwd : instance of Forward - The forward solution used to create the source estimate. - cov : instance of Covariance - The noise covariance used to estimate the resting cortical - activations. Should be an evoked covariance, not empty room. - - Returns - ------- - snr_stc : instance of SourceEstimate - The source estimate with the SNR computed. - - Notes - ----- - We define the SNR in decibels for each source location at each - time point as: - - .. math:: - - {\\rm SNR} = 10\\log_10[\\frac{a^2}{N}\\sum_k\\frac{b_k^2}{s_k^2}] - - where :math:`\\b_k` is the signal on sensor :math:`k` provided by the - forward model for a source with unit amplitude, :math:`a` is the - source amplitude, :math:`N` is the number of sensors, and - :math:`s_k^2` is the noise variance on sensor :math:`k`. - """ - from .forward import convert_forward_solution - from .minimum_norm.inverse import _prepare_forward - - fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - print("done convert forward solution\n") - - - - #def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, - # use_cps, exp, limit_depth_chs, combine_xyz, - # allow_fixed_depth, limit): - -# return (forward, info_picked, gain, depth_prior, orient_prior, source_std, -# trace_GRGT, noise_cov, whitener) - - _, _, _, _, _, A, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, - use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', - allow_fixed_depth=True, limit=None) - - print("done _prepare_forward\n") - N = cov['dim'] - print ("N = %d sensors\n" % N) - b_k2 = (A * A).T - s_k2 = np.diag(cov['data']) - - print (np.shape(b_k2)) - print (np.shape(s_k2)) - scaling = (1. / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] -# print (np.shape(scaling)) - - print("End of SNR Func\n") - \ No newline at end of file diff --git a/tutorials/source-modeling/plot_mne_solutions.py b/tutorials/source-modeling/plot_mne_solutions.py index 8050e3fc98f..3a9faf24c6f 100644 --- a/tutorials/source-modeling/plot_mne_solutions.py +++ b/tutorials/source-modeling/plot_mne_solutions.py @@ -49,6 +49,10 @@ size=(600, 600)) stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) +stc.estimate_snr(evoked.info, fwd, cov) + +sdfsdfdsf + brain = stc.plot(figure=1, **kwargs) brain.add_text(0.1, 0.9, 'MNE', 'title', font_size=14) From f9fb26bd568ac151079e05a3d3d50740e6db2d52 Mon Sep 17 00:00:00 2001 From: Mainak Jas Date: Tue, 1 Oct 2019 15:27:09 -0400 Subject: [PATCH 03/24] ENH: name change --- mne/source_estimate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index bc6927f4149..94c333129ea 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1157,14 +1157,14 @@ def estimate_snr(self, info, fwd, cov, verbose=None): # return (forward, info_picked, gain, depth_prior, orient_prior, source_std, # trace_GRGT, noise_cov, whitener) - _, _, _, _, _, A, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, + _, _, _, _, _, source_std, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', allow_fixed_depth=True, limit=None) print("done _prepare_forward\n") N = cov['dim'] print ("N = %d sensors\n" % N) - b_k2 = (A * A).T + b_k2 = (source_std * source_std).T s_k2 = np.diag(cov['data']) print (np.shape(b_k2)) From d954ae4915b17b45e45f4f5daaef95642afe1937 Mon Sep 17 00:00:00 2001 From: Mainak Jas Date: Tue, 1 Oct 2019 15:34:46 -0400 Subject: [PATCH 04/24] Make it work --- mne/source_estimate.py | 9 ++++++--- tutorials/source-modeling/plot_mne_solutions.py | 4 +++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 94c333129ea..197737b50df 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1145,6 +1145,7 @@ def estimate_snr(self, info, fwd, cov, verbose=None): from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward + snr_stc = self.copy() fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) print("done convert forward solution\n") @@ -1157,22 +1158,24 @@ def estimate_snr(self, info, fwd, cov, verbose=None): # return (forward, info_picked, gain, depth_prior, orient_prior, source_std, # trace_GRGT, noise_cov, whitener) - _, _, _, _, _, source_std, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, + _, _, G, _, _, source_std, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', allow_fixed_depth=True, limit=None) print("done _prepare_forward\n") N = cov['dim'] print ("N = %d sensors\n" % N) - b_k2 = (source_std * source_std).T + b_k2 = (G * G).T s_k2 = np.diag(cov['data']) print (np.shape(b_k2)) print (np.shape(s_k2)) scaling = (1. / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] # print (np.shape(scaling)) - + snr_stc._data = (snr_stc.data * snr_stc.data) * scaling + print("End of SNR Func\n") + return snr_stc def _center_of_mass(vertices, values, hemi, surf, subject, subjects_dir, diff --git a/tutorials/source-modeling/plot_mne_solutions.py b/tutorials/source-modeling/plot_mne_solutions.py index 3a9faf24c6f..fd512b92297 100644 --- a/tutorials/source-modeling/plot_mne_solutions.py +++ b/tutorials/source-modeling/plot_mne_solutions.py @@ -49,7 +49,9 @@ size=(600, 600)) stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) -stc.estimate_snr(evoked.info, fwd, cov) +stc_snr = stc.estimate_snr(evoked.info, fwd, cov) + +stc_snr.plot(subjects_dir=subjects_dir) sdfsdfdsf From aef9f00a43e2be3737d7d100560b4b7a356c76fe Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Tue, 1 Oct 2019 17:22:00 -0400 Subject: [PATCH 05/24] SNR works, now clean up code --- mne/source_estimate.py | 31 ++++++++----------- .../source-modeling/plot_mne_solutions.py | 8 +++-- 2 files changed, 18 insertions(+), 21 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 197737b50df..121c177b59c 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1147,13 +1147,6 @@ def estimate_snr(self, info, fwd, cov, verbose=None): snr_stc = self.copy() fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - print("done convert forward solution\n") - - - - #def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, - # use_cps, exp, limit_depth_chs, combine_xyz, - # allow_fixed_depth, limit): # return (forward, info_picked, gain, depth_prior, orient_prior, source_std, # trace_GRGT, noise_cov, whitener) @@ -1161,19 +1154,21 @@ def estimate_snr(self, info, fwd, cov, verbose=None): _, _, G, _, _, source_std, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', allow_fixed_depth=True, limit=None) - - print("done _prepare_forward\n") - N = cov['dim'] - print ("N = %d sensors\n" % N) - b_k2 = (G * G).T + nTimes = np.shape(snr_stc.data)[1] + + N = cov['dim'] # number of sensors/channels + # G is gain matrix [channels x sources] + M = np.shape(snr_stc.data)[0] # number of sources s_k2 = np.diag(cov['data']) - print (np.shape(b_k2)) - print (np.shape(s_k2)) - scaling = (1. / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] - # print (np.shape(scaling)) - snr_stc._data = (snr_stc.data * snr_stc.data) * scaling - + for t in np.arange(0,nTimes): # for each time pt of the data + print ('t = %d\n' % t) + for i in np.arange(0,M): # loop over sources + b_k2 = G[:,i]*G[:,i] #ith source + scaling = (1 / N) * np.sum(b_k2 / s_k2) # scalar #, axis=1)[:, np.newaxis] + a_it = snr_stc.data[i,t]*snr_stc.data[i,t] + snr_stc.data[i,t] = 10*np.log10(a_it * scaling) + print("End of SNR Func\n") return snr_stc diff --git a/tutorials/source-modeling/plot_mne_solutions.py b/tutorials/source-modeling/plot_mne_solutions.py index fd512b92297..eff120a62f9 100644 --- a/tutorials/source-modeling/plot_mne_solutions.py +++ b/tutorials/source-modeling/plot_mne_solutions.py @@ -49,13 +49,15 @@ size=(600, 600)) stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) -stc_snr = stc.estimate_snr(evoked.info, fwd, cov) +snr_stc = stc.estimate_snr(evoked.info, fwd, cov) -stc_snr.plot(subjects_dir=subjects_dir) +kwargs2 = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, + size=(600, 600),colormap='coolwarm',clim=dict(kind='value', lims=(-20,0,20))) +snr_stc.plot(figure=1, **kwargs2) sdfsdfdsf -brain = stc.plot(figure=1, **kwargs) +brain = stc.plot(figure=3, **kwargs) brain.add_text(0.1, 0.9, 'MNE', 'title', font_size=14) From a013fe784cda60934994d6834f9700a98f844276 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 10:29:08 -0400 Subject: [PATCH 06/24] fixed unnecessary stc copy, returns stc object --- mne/source_estimate.py | 52 +++++++++++++++++++++--------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 121c177b59c..2d9852299b5 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -25,8 +25,6 @@ plot_volume_source_estimates) from .io.base import ToDataFrameMixin, TimeMixin from .externals.h5io import read_hdf5, write_hdf5 -#from mne.forward import convert_forward_solution -#from mne.minimum_norm.inverse import _prepare_forward def _read_stc(filename): """Aux Function.""" @@ -1101,7 +1099,7 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False): @verbose def estimate_snr(self, info, fwd, cov, verbose=None): - """Compute time-varying SNR + """Compute time-varying SNR in the source space This function should only be used with source estimates with units nA (i.e., MNE-like solutions, *not* dSPM or sLORETA). @@ -1144,32 +1142,34 @@ def estimate_snr(self, info, fwd, cov, verbose=None): """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward - - snr_stc = self.copy() - fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - - # return (forward, info_picked, gain, depth_prior, orient_prior, source_std, - # trace_GRGT, noise_cov, whitener) + from mne import SourceEstimate - _, _, G, _, _, source_std, _, cov, _ = _prepare_forward(fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, - use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', - allow_fixed_depth=True, limit=None) - nTimes = np.shape(snr_stc.data)[1] + M = np.shape(self.data)[0] + T = np.shape(self.data)[1] + snr_stc = SourceEstimate(np.zeros((M, T)), self.vertices, self.tmin, + self.tstep, self.subject) - N = cov['dim'] # number of sensors/channels - # G is gain matrix [channels x sources] - M = np.shape(snr_stc.data)[0] # number of sources - s_k2 = np.diag(cov['data']) + fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - for t in np.arange(0,nTimes): # for each time pt of the data - print ('t = %d\n' % t) - for i in np.arange(0,M): # loop over sources - b_k2 = G[:,i]*G[:,i] #ith source - scaling = (1 / N) * np.sum(b_k2 / s_k2) # scalar #, axis=1)[:, np.newaxis] - a_it = snr_stc.data[i,t]*snr_stc.data[i,t] - snr_stc.data[i,t] = 10*np.log10(a_it * scaling) - - print("End of SNR Func\n") + # G is gain matrix [ch x src], cov is noise covariance [ch x ch] + _, _, G, _, _, _, _, cov, _ = _prepare_forward(fwd, info, cov, + fixed=True, + loose=0, + rank=None, + pca=False, + use_cps=True, + exp=None, + limit_depth_chs=False, + combine_xyz='fro', + allow_fixed_depth=True, + limit=None) + + N = cov['dim'] # number of sensors/channels + b_k2 = (G * G).T + s_k2 = np.diag(cov['data']) + scaling = (1/N) * np.sum(b_k2 /s_k2, axis=1)[:, np.newaxis] + snr_stc.data = 10*np.log10((self.data * self.data) * scaling) + return snr_stc From 3b4a93d720356c77d8659bd1e6f9c10655fcfd14 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 10:46:38 -0400 Subject: [PATCH 07/24] added plot source space SNR tutorial file --- .../source-modeling/plot_source_space_SNR.py | 56 +++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 tutorials/source-modeling/plot_source_space_SNR.py diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py new file mode 100644 index 00000000000..1b03b2928e6 --- /dev/null +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -0,0 +1,56 @@ +# -*- coding: utf-8 -*- +""" +.. _tut-mne-fixed-free: + +=============================== +Computing source space SNR +=============================== + +This example shows example fixed- and free-orientation source localizations +produced by MNE, dSPM, sLORETA, and eLORETA. +""" +# Author: Eric Larson +# +# License: BSD (3-clause) + +import mne +from mne.datasets import sample +from mne.minimum_norm import make_inverse_operator, apply_inverse + +print(__doc__) + +data_path = sample.data_path() +subjects_dir = data_path + '/subjects' + +# Read data +fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif' +evoked = mne.read_evokeds(fname_evoked, condition='Left Auditory', + baseline=(None, 0)) +fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif' +fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif' +fwd = mne.read_forward_solution(fname_fwd) +cov = mne.read_cov(fname_cov) + +############################################################################### +# Fixed orientation +# ----------------- +# First let's create a fixed-orientation inverse, with the default weighting. + +inv = make_inverse_operator(evoked.info, fwd, cov, loose=0., depth=0.8, + verbose=True) + +############################################################################### +# Let's look at the current estimates using MNE. We'll take the absolute +# value of the source estimates to simplify the visualization. + +snr = 3.0 +lambda2 = 1.0 / snr ** 2 +kwargs = dict(initial_time=0.08, hemi='both', subjects_dir=subjects_dir, + size=(600, 600)) + +stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) +snr_stc = stc.estimate_snr(evoked.info, fwd, cov) + +kwargs2 = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, + size=(600, 600),colormap='coolwarm',clim=dict(kind='value', lims=(-20,0,20))) +snr_stc.plot(figure=1, **kwargs2) From a80c567fc20a72ea88f14af40507fa9091aaa9de Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 10:57:29 -0400 Subject: [PATCH 08/24] updated author names --- tutorials/source-modeling/plot_mne_solutions.py | 10 +--------- tutorials/source-modeling/plot_source_space_SNR.py | 3 ++- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/tutorials/source-modeling/plot_mne_solutions.py b/tutorials/source-modeling/plot_mne_solutions.py index eff120a62f9..8050e3fc98f 100644 --- a/tutorials/source-modeling/plot_mne_solutions.py +++ b/tutorials/source-modeling/plot_mne_solutions.py @@ -49,15 +49,7 @@ size=(600, 600)) stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) -snr_stc = stc.estimate_snr(evoked.info, fwd, cov) - -kwargs2 = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, - size=(600, 600),colormap='coolwarm',clim=dict(kind='value', lims=(-20,0,20))) -snr_stc.plot(figure=1, **kwargs2) - -sdfsdfdsf - -brain = stc.plot(figure=3, **kwargs) +brain = stc.plot(figure=1, **kwargs) brain.add_text(0.1, 0.9, 'MNE', 'title', font_size=14) diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index 1b03b2928e6..98645c593fa 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -9,7 +9,8 @@ This example shows example fixed- and free-orientation source localizations produced by MNE, dSPM, sLORETA, and eLORETA. """ -# Author: Eric Larson +# Author: Padma Sundaram +# Kaisu Lankinen # # License: BSD (3-clause) From 5ad490f02ee36fde0547bf4d1e0754fe0122f006 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:12:43 -0400 Subject: [PATCH 09/24] fixed period for pydocsty and suggested documenation edit --- mne/source_estimate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 2d9852299b5..43456cffe7d 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1099,10 +1099,10 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False): @verbose def estimate_snr(self, info, fwd, cov, verbose=None): - """Compute time-varying SNR in the source space + """Compute time-varying SNR in the source space. This function should only be used with source estimates with units - nA (i.e., MNE-like solutions, *not* dSPM or sLORETA). + nanoAmperes (i.e., MNE-like solutions, *not* dSPM or sLORETA). Reference: Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., From 9bf444e273b31cd97f7a8d5918dfde70263d23f3 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:26:16 -0400 Subject: [PATCH 10/24] edits suggested from @drammock --- mne/source_estimate.py | 32 ++++++++----------- .../source-modeling/plot_source_space_SNR.py | 5 ++- 2 files changed, 15 insertions(+), 22 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 43456cffe7d..e206de5dcda 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1099,18 +1099,11 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False): @verbose def estimate_snr(self, info, fwd, cov, verbose=None): - """Compute time-varying SNR in the source space. + r"""Compute time-varying SNR in the source space. This function should only be used with source estimates with units nanoAmperes (i.e., MNE-like solutions, *not* dSPM or sLORETA). - Reference: - Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., - Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). - Mapping the Signal-To-Noise-Ratios of Cortical Sources in - Magnetoencephalography and Electroencephalography. - Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 - Parameters ---------- info : instance of io.meas_info.Info @@ -1133,12 +1126,19 @@ def estimate_snr(self, info, fwd, cov, verbose=None): .. math:: - {\\rm SNR} = 10\\log_10[\\frac{a^2}{N}\\sum_k\\frac{b_k^2}{s_k^2}] + {\rm SNR} = 10\log_10[\frac{a^2}{N}\sum_k\frac{b_k^2}{s_k^2}] where :math:`\\b_k` is the signal on sensor :math:`k` provided by the forward model for a source with unit amplitude, :math:`a` is the source amplitude, :math:`N` is the number of sensors, and :math:`s_k^2` is the noise variance on sensor :math:`k`. + + Reference: + Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., + Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + Mapping the Signal-To-Noise-Ratios of Cortical Sources in + Magnetoencephalography and Electroencephalography. + Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward @@ -1153,16 +1153,10 @@ def estimate_snr(self, info, fwd, cov, verbose=None): # G is gain matrix [ch x src], cov is noise covariance [ch x ch] _, _, G, _, _, _, _, cov, _ = _prepare_forward(fwd, info, cov, - fixed=True, - loose=0, - rank=None, - pca=False, - use_cps=True, - exp=None, - limit_depth_chs=False, - combine_xyz='fro', - allow_fixed_depth=True, - limit=None) + fixed=True, loose=0, rank=None, pca=False, + use_cps=True, exp=None, limit_depth_chs=False, + combine_xyz='fro', allow_fixed_depth=True, + limit=None) N = cov['dim'] # number of sensors/channels b_k2 = (G * G).T diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index 98645c593fa..171da8544ad 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -1,13 +1,12 @@ # -*- coding: utf-8 -*- """ -.. _tut-mne-fixed-free: +.. _tut-mne-sourcespace-snr: =============================== Computing source space SNR =============================== -This example shows example fixed- and free-orientation source localizations -produced by MNE, dSPM, sLORETA, and eLORETA. +This example shows how to compute and plot source space SNR. """ # Author: Padma Sundaram # Kaisu Lankinen From 6e3106b0e02a45447e4013bc4d65f4ac127dfb10 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:47:46 -0400 Subject: [PATCH 11/24] made edits suggested by reviewers --- mne/source_estimate.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index e206de5dcda..d1d4cfb6a4c 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1142,7 +1142,6 @@ def estimate_snr(self, info, fwd, cov, verbose=None): """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward - from mne import SourceEstimate M = np.shape(self.data)[0] T = np.shape(self.data)[1] @@ -1158,11 +1157,11 @@ def estimate_snr(self, info, fwd, cov, verbose=None): combine_xyz='fro', allow_fixed_depth=True, limit=None) - N = cov['dim'] # number of sensors/channels + N = info.cov['dim'] # number of sensors/channels b_k2 = (G * G).T s_k2 = np.diag(cov['data']) - scaling = (1/N) * np.sum(b_k2 /s_k2, axis=1)[:, np.newaxis] - snr_stc.data = 10*np.log10((self.data * self.data) * scaling) + scaling = (1 / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] + snr_stc.data = 10 * np.log10((self.data * self.data) * scaling) return snr_stc From 04012cfb19c96a2f708ac2cb080cfbada469faf4 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:49:21 -0400 Subject: [PATCH 12/24] made edits suggested by reviewers --- mne/source_estimate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index d1d4cfb6a4c..400ed0883c6 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1157,7 +1157,7 @@ def estimate_snr(self, info, fwd, cov, verbose=None): combine_xyz='fro', allow_fixed_depth=True, limit=None) - N = info.cov['dim'] # number of sensors/channels + N = cov['dim'] # number of sensors/channels b_k2 = (G * G).T s_k2 = np.diag(cov['data']) scaling = (1 / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] From 3213d2bd052a4c418c8e25493d9ef3956c830c29 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:51:25 -0400 Subject: [PATCH 13/24] N -> n_channels --- mne/source_estimate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 400ed0883c6..55d150a7fab 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1157,10 +1157,10 @@ def estimate_snr(self, info, fwd, cov, verbose=None): combine_xyz='fro', allow_fixed_depth=True, limit=None) - N = cov['dim'] # number of sensors/channels + n_channels = cov['dim'] # number of sensors/channels b_k2 = (G * G).T s_k2 = np.diag(cov['data']) - scaling = (1 / N) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] + scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] snr_stc.data = 10 * np.log10((self.data * self.data) * scaling) return snr_stc From 3efeba957c74f98286c58634e6fefdc9d2c8e076 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:54:50 -0400 Subject: [PATCH 14/24] Fixed formatting of Reference in comment --- mne/source_estimate.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 55d150a7fab..adbaaa07d3b 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1132,10 +1132,11 @@ def estimate_snr(self, info, fwd, cov, verbose=None): forward model for a source with unit amplitude, :math:`a` is the source amplitude, :math:`N` is the number of sensors, and :math:`s_k^2` is the noise variance on sensor :math:`k`. - - Reference: - Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., - Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + + References: + ---------- + .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, + D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). Mapping the Signal-To-Noise-Ratios of Cortical Sources in Magnetoencephalography and Electroencephalography. Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 From beac77246f87a73647453a743f393bee227ea7f4 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 11:58:07 -0400 Subject: [PATCH 15/24] fixed formatting in doc and in code --- mne/source_estimate.py | 7 +++---- tutorials/source-modeling/plot_source_space_SNR.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index adbaaa07d3b..a0588860d5e 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1143,10 +1143,9 @@ def estimate_snr(self, info, fwd, cov, verbose=None): """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward - - M = np.shape(self.data)[0] - T = np.shape(self.data)[1] - snr_stc = SourceEstimate(np.zeros((M, T)), self.vertices, self.tmin, + + n_src, n_times = self.data.shape + snr_stc = SourceEstimate(np.zeros((n_src, n_times)), self.vertices, self.tmin, self.tstep, self.subject) fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index 171da8544ad..a1438554137 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -9,7 +9,7 @@ This example shows how to compute and plot source space SNR. """ # Author: Padma Sundaram -# Kaisu Lankinen +# Kaisu Lankinen # # License: BSD (3-clause) From 60c1fd0a0a294ef8556dc2d832b86295947f0cbc Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 12:16:17 -0400 Subject: [PATCH 16/24] fixed snr_stc creation --- mne/source_estimate.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index a0588860d5e..15ee8948948 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1143,10 +1143,6 @@ def estimate_snr(self, info, fwd, cov, verbose=None): """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward - - n_src, n_times = self.data.shape - snr_stc = SourceEstimate(np.zeros((n_src, n_times)), self.vertices, self.tmin, - self.tstep, self.subject) fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) @@ -1161,7 +1157,8 @@ def estimate_snr(self, info, fwd, cov, verbose=None): b_k2 = (G * G).T s_k2 = np.diag(cov['data']) scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] - snr_stc.data = 10 * np.log10((self.data * self.data) * scaling) + snr_stc = SourceEstimate(10 * np.log10((self.data * self.data) * scaling), + self.vertices, self.tmin, self.tstep, self.subject) return snr_stc From dbcff92318bf0f79502d3e176bd349596baf4d1b Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Wed, 2 Oct 2019 14:58:24 -0400 Subject: [PATCH 17/24] further minor edits --- mne/source_estimate.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 15ee8948948..d0eaab90247 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -23,8 +23,8 @@ fill_doc, _check_option, _validate_type, _check_src_normal) from .viz import (plot_source_estimates, plot_vector_source_estimates, plot_volume_source_estimates) -from .io.base import ToDataFrameMixin, TimeMixin -from .externals.h5io import read_hdf5, write_hdf5 +from .io.base import (ToDataFrameMixin, TimeMixin) +from .externals.h5io import (read_hdf5, write_hdf5) def _read_stc(filename): """Aux Function.""" @@ -1136,10 +1136,10 @@ def estimate_snr(self, info, fwd, cov, verbose=None): References: ---------- .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, - D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). - Mapping the Signal-To-Noise-Ratios of Cortical Sources in - Magnetoencephalography and Electroencephalography. - Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 + D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + Mapping the Signal-To-Noise-Ratios of Cortical Sources in + Magnetoencephalography and Electroencephalography. + Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward From 66fdd89bf527298735f5d9d6f174a61b4f4c664c Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Thu, 3 Oct 2019 09:09:47 -0400 Subject: [PATCH 18/24] fixed formatting issues and colormap choice in eg --- mne/source_estimate.py | 5 +++-- tutorials/source-modeling/plot_source_space_SNR.py | 11 +++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index d0eaab90247..26ba9a0b79a 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -23,8 +23,9 @@ fill_doc, _check_option, _validate_type, _check_src_normal) from .viz import (plot_source_estimates, plot_vector_source_estimates, plot_volume_source_estimates) -from .io.base import (ToDataFrameMixin, TimeMixin) -from .externals.h5io import (read_hdf5, write_hdf5) +from .io.base import ToDataFrameMixin, TimeMixin +from .externals.h5io import read_hdf5, write_hdf5 + def _read_stc(filename): """Aux Function.""" diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index a1438554137..16f9381e5d1 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -45,12 +45,11 @@ snr = 3.0 lambda2 = 1.0 / snr ** 2 -kwargs = dict(initial_time=0.08, hemi='both', subjects_dir=subjects_dir, - size=(600, 600)) -stc = abs(apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True)) +stc = apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True) snr_stc = stc.estimate_snr(evoked.info, fwd, cov) -kwargs2 = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, - size=(600, 600),colormap='coolwarm',clim=dict(kind='value', lims=(-20,0,20))) -snr_stc.plot(figure=1, **kwargs2) +kwargs = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, + size=(600, 600), colormap='viridis', + clim=dict(kind='value', lims=(-20, 0, 20))) +snr_stc.plot(figure=1, **kwargs) From 39b64fba73d6c5ed2fb76aed272590eb1657ab72 Mon Sep 17 00:00:00 2001 From: Padma Sundaram Date: Thu, 3 Oct 2019 11:14:21 -0400 Subject: [PATCH 19/24] removed colon after Reference --- mne/source_estimate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 26ba9a0b79a..0d4e5086234 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1134,7 +1134,7 @@ def estimate_snr(self, info, fwd, cov, verbose=None): source amplitude, :math:`N` is the number of sensors, and :math:`s_k^2` is the noise variance on sensor :math:`k`. - References: + References ---------- .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). From 48a7c706e4b1ab4ea754ffb635aae47ec1738066 Mon Sep 17 00:00:00 2001 From: Kaisu Lankinen Date: Thu, 3 Oct 2019 16:44:37 -0400 Subject: [PATCH 20/24] Corrected SNR calculation to take unwhitened sensor signal b instead of whitened, updated tutorial file --- mne/source_estimate.py | 141 ++++++++++-------- .../source-modeling/plot_source_space_SNR.py | 103 +++++++++++-- 2 files changed, 165 insertions(+), 79 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 0d4e5086234..6b2975f6fb9 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1098,71 +1098,6 @@ def transform(self, func, idx=None, tmin=None, tmax=None, copy=False): return stcs - @verbose - def estimate_snr(self, info, fwd, cov, verbose=None): - r"""Compute time-varying SNR in the source space. - - This function should only be used with source estimates with units - nanoAmperes (i.e., MNE-like solutions, *not* dSPM or sLORETA). - - Parameters - ---------- - info : instance of io.meas_info.Info - The measurement info. - fwd : instance of Forward - The forward solution used to create the source estimate. - cov : instance of Covariance - The noise covariance used to estimate the resting cortical - activations. Should be an evoked covariance, not empty room. - - Returns - ------- - snr_stc : instance of SourceEstimate - The source estimate with the SNR computed. - - Notes - ----- - We define the SNR in decibels for each source location at each - time point as: - - .. math:: - - {\rm SNR} = 10\log_10[\frac{a^2}{N}\sum_k\frac{b_k^2}{s_k^2}] - - where :math:`\\b_k` is the signal on sensor :math:`k` provided by the - forward model for a source with unit amplitude, :math:`a` is the - source amplitude, :math:`N` is the number of sensors, and - :math:`s_k^2` is the noise variance on sensor :math:`k`. - - References - ---------- - .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, - D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). - Mapping the Signal-To-Noise-Ratios of Cortical Sources in - Magnetoencephalography and Electroencephalography. - Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 - """ - from .forward import convert_forward_solution - from .minimum_norm.inverse import _prepare_forward - - fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - - # G is gain matrix [ch x src], cov is noise covariance [ch x ch] - _, _, G, _, _, _, _, cov, _ = _prepare_forward(fwd, info, cov, - fixed=True, loose=0, rank=None, pca=False, - use_cps=True, exp=None, limit_depth_chs=False, - combine_xyz='fro', allow_fixed_depth=True, - limit=None) - - n_channels = cov['dim'] # number of sensors/channels - b_k2 = (G * G).T - s_k2 = np.diag(cov['data']) - scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1)[:, np.newaxis] - snr_stc = SourceEstimate(10 * np.log10((self.data * self.data) * scaling), - self.vertices, self.tmin, self.tstep, self.subject) - - return snr_stc - def _center_of_mass(vertices, values, hemi, surf, subject, subjects_dir, restrict_vertices): @@ -1570,6 +1505,81 @@ def extract_label_time_course(self, labels, src, mode='mean_flip', return label_tc + + @verbose + def estimate_snr(self, info, fwd, cov, verbose=None): + r"""Compute time-varying SNR in the source space. + + This function should only be used with source estimates with units + nanoAmperes (i.e., MNE-like solutions, *not* dSPM or sLORETA). + + .. warning:: This function currently only works properly for fixed + orientation... + + Parameters + ---------- + info : instance of io.meas_info.Info + The measurement info. + fwd : instance of Forward + The forward solution used to create the source estimate. + cov : instance of Covariance + The noise covariance used to estimate the resting cortical + activations. Should be an evoked covariance, not empty room. + + Returns + ------- + snr_stc : instance of SourceEstimate + The source estimate with the SNR computed. + + Notes + ----- + We define the SNR in decibels for each source location at each + time point as: + + .. math:: + + {\rm SNR} = 10\log_10[\frac{a^2}{N}\sum_k\frac{b_k^2}{s_k^2}] + + where :math:`\\b_k` is the signal on sensor :math:`k` provided by the + forward model for a source with unit amplitude, :math:`a` is the + source amplitude, :math:`N` is the number of sensors, and + :math:`s_k^2` is the noise variance on sensor :math:`k`. + + References + ---------- + .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, + D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). + Mapping the Signal-To-Noise-Ratios of Cortical Sources in + Magnetoencephalography and Electroencephalography. + Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 + """ + from .forward import convert_forward_solution + from .minimum_norm.inverse import _prepare_forward + + if (self.data >= 0).all(): + warn_('This STC appears to be from free orientation, currently SNR ' + 'function is valid only for fixed orientation') + + fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) + + # G is gain matrix [ch x src], cov is noise covariance [ch x ch] + G, _, _, _, _, _, _, cov, _ = _prepare_forward(fwd, info, cov, + fixed=True, loose=0, rank=None, pca=False, + use_cps=True, exp=None, limit_depth_chs=False, + combine_xyz='fro', allow_fixed_depth=True, + limit=None) + G = G['sol']['data'] + n_channels = cov['dim'] # number of sensors/channels + print('Number of channels in calculating SNR: %d' % n_channels) + b_k2 = (G * G).T + s_k2 = np.diag(cov['data']) + scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1, keepdims=True) + snr_stc = SourceEstimate(10 * np.log10((self.data * self.data) * scaling), + self.vertices, self.tmin, self.tstep, self.subject) + + return snr_stc + + def get_peak(self, hemi=None, tmin=None, tmax=None, mode='abs', vert_as_index=False, time_as_index=False): """Get location and latency of peak amplitude. @@ -1703,6 +1713,7 @@ def center_of_mass(self, subject=None, hemi=None, restrict_vertices=False, return vertex, hemi, t + class _BaseVectorSourceEstimate(_BaseSourceEstimate): _data_ndim = 3 diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index 16f9381e5d1..2e0496b5474 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -7,6 +7,13 @@ =============================== This example shows how to compute and plot source space SNR. + +Reference: +Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., Ishitobi, M., +Vaina, L. M., & Stufflebeam, S. M. (2009). Mapping the Signal-To-Noise-Ratios +of Cortical Sources in Magnetoencephalography and Electroencephalography. +Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 + """ # Author: Padma Sundaram # Kaisu Lankinen @@ -15,7 +22,9 @@ import mne from mne.datasets import sample -from mne.minimum_norm import make_inverse_operator, apply_inverse +from mne.minimum_norm import make_inverse_operator, apply_inverse, read_inverse_operator +import numpy as np +import matplotlib.pyplot as plt print(__doc__) @@ -32,24 +41,90 @@ cov = mne.read_cov(fname_cov) ############################################################################### -# Fixed orientation -# ----------------- -# First let's create a fixed-orientation inverse, with the default weighting. - -inv = make_inverse_operator(evoked.info, fwd, cov, loose=0., depth=0.8, - verbose=True) +# MEG-EEG: -############################################################################### -# Let's look at the current estimates using MNE. We'll take the absolute -# value of the source estimates to simplify the visualization. +# Read inverse operator: +inv_op = make_inverse_operator(evoked.info, fwd, cov, fixed=True, verbose=True) +# Calculate MNE: snr = 3.0 lambda2 = 1.0 / snr ** 2 +stc = apply_inverse(evoked, inv_op, lambda2, 'MNE', verbose=True) -stc = apply_inverse(evoked, inv, lambda2, 'MNE', verbose=True) +# Calculate SNR in source space: snr_stc = stc.estimate_snr(evoked.info, fwd, cov) -kwargs = dict(initial_time=0.29, hemi='both', subjects_dir=subjects_dir, - size=(600, 600), colormap='viridis', - clim=dict(kind='value', lims=(-20, 0, 20))) +# Plot an average SNR across source points over time: +ave = np.mean(snr_stc.data, axis=0) + +fig, ax = plt.subplots() +ax.plot(evoked.times, ave) +ax.set(xlabel='Time (sec)', ylabel='SNR MEG-EEG') +fig.tight_layout() + +# Find time point of maximum SNR: +maxidx = np.argmax(ave) + +# Plot SNR on source space at the time point of maximum SNR: +kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', + subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', + clim=dict(kind='value', lims=(-100, -75, -50))) snr_stc.plot(figure=1, **kwargs) + +############################################################################### +# EEG: + +evoked_eeg = evoked.copy().pick_types(eeg=True, meg=False) + +# Make inverse oparator: +inv_op_eeg = make_inverse_operator(evoked_eeg.info, fwd, cov, fixed=True, + verbose=True) + +# Calculate MNE: +stc_eeg = apply_inverse(evoked_eeg, inv_op_eeg, lambda2, 'MNE', verbose=True) + +# Calculate SNR in source space: +snr_stc_eeg = stc_eeg.estimate_snr(evoked_eeg.info, fwd, cov) + +# Plot an average SNR across source points over time: +ave_eeg = np.mean(snr_stc_eeg.data, axis=0) +fig, ax = plt.subplots() +ax.plot(evoked_eeg.times, ave_eeg) +ax.set(xlabel='Time (sec)', ylabel='SNR EEG') +fig.tight_layout() + +# Plot SNR on source space at the time point of maximum SNR (same time point as +# for MEG-EEG) +kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', + subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', + clim=dict(kind='value', lims=(-100, -75, -50))) +snr_stc_eeg.plot(figure=2, **kwargs) + +############################################################################### +# MEG: + +evoked_meg = evoked.copy().pick_types(eeg=False, meg=True) + +# Make inverse oparator: +inv_op_meg = make_inverse_operator(evoked_meg.info, fwd, cov, fixed=True, + verbose=True) + +# Calculate MNE: +stc_meg = apply_inverse(evoked_meg, inv_op_meg, lambda2, 'MNE', verbose=True) + +# Calculate SNR in source space: +snr_stc_meg = stc_meg.estimate_snr(evoked_meg.info, fwd, cov) + +# Plot an average SNR across source points over time: +ave_meg = np.mean(snr_stc_meg.data, axis=0) +fig, ax = plt.subplots() +ax.plot(evoked_meg.times, ave_meg) +ax.set(xlabel='Time (sec)', ylabel='SNR MEG') +fig.tight_layout() + +# Plot SNR on source space at the time point of maximum SNR (same time point as +# for MEG-EEG) +kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', + subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', + clim=dict(kind='value', lims=(-100, -75, -50))) +snr_stc_meg.plot(figure=3, **kwargs) From 7cd00576f42fb890b4667bd4280aef04a957796d Mon Sep 17 00:00:00 2001 From: Kaisu Lankinen Date: Fri, 4 Oct 2019 10:56:52 -0400 Subject: [PATCH 21/24] Minor changes to formatting --- mne/source_estimate.py | 5 +++-- tutorials/source-modeling/plot_source_space_SNR.py | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 6b2975f6fb9..43025f0d0d3 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1514,17 +1514,18 @@ def estimate_snr(self, info, fwd, cov, verbose=None): nanoAmperes (i.e., MNE-like solutions, *not* dSPM or sLORETA). .. warning:: This function currently only works properly for fixed - orientation... + orientation. Parameters ---------- - info : instance of io.meas_info.Info + info : instance Info The measurement info. fwd : instance of Forward The forward solution used to create the source estimate. cov : instance of Covariance The noise covariance used to estimate the resting cortical activations. Should be an evoked covariance, not empty room. + %(verbose)s Returns ------- diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/tutorials/source-modeling/plot_source_space_SNR.py index 2e0496b5474..2a5dcc31ed8 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/tutorials/source-modeling/plot_source_space_SNR.py @@ -20,6 +20,8 @@ # # License: BSD (3-clause) +# sphinx_gallery_thumbnail_number = 2 + import mne from mne.datasets import sample from mne.minimum_norm import make_inverse_operator, apply_inverse, read_inverse_operator From fb2655fc02c066b028ca3b03849d25ab78f9baf2 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 4 Oct 2019 12:54:05 -0400 Subject: [PATCH 22/24] FIX: transparent --- .../inverse}/plot_source_space_SNR.py | 84 ++++++------------- mne/source_estimate.py | 39 ++++----- 2 files changed, 44 insertions(+), 79 deletions(-) rename {tutorials/source-modeling => examples/inverse}/plot_source_space_SNR.py (55%) diff --git a/tutorials/source-modeling/plot_source_space_SNR.py b/examples/inverse/plot_source_space_SNR.py similarity index 55% rename from tutorials/source-modeling/plot_source_space_SNR.py rename to examples/inverse/plot_source_space_SNR.py index 2a5dcc31ed8..9813aca28ef 100644 --- a/tutorials/source-modeling/plot_source_space_SNR.py +++ b/examples/inverse/plot_source_space_SNR.py @@ -1,19 +1,12 @@ # -*- coding: utf-8 -*- """ -.. _tut-mne-sourcespace-snr: +.. ex-mne-sourcespace-snr: =============================== Computing source space SNR =============================== -This example shows how to compute and plot source space SNR. - -Reference: -Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., Ishitobi, M., -Vaina, L. M., & Stufflebeam, S. M. (2009). Mapping the Signal-To-Noise-Ratios -of Cortical Sources in Magnetoencephalography and Electroencephalography. -Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 - +This example shows how to compute and plot source space SNR as in [1]_. """ # Author: Padma Sundaram # Kaisu Lankinen @@ -24,7 +17,7 @@ import mne from mne.datasets import sample -from mne.minimum_norm import make_inverse_operator, apply_inverse, read_inverse_operator +from mne.minimum_norm import make_inverse_operator, apply_inverse import numpy as np import matplotlib.pyplot as plt @@ -43,7 +36,8 @@ cov = mne.read_cov(fname_cov) ############################################################################### -# MEG-EEG: +# MEG-EEG +# ------- # Read inverse operator: inv_op = make_inverse_operator(evoked.info, fwd, cov, fixed=True, verbose=True) @@ -68,65 +62,41 @@ maxidx = np.argmax(ave) # Plot SNR on source space at the time point of maximum SNR: -kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', - subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', - clim=dict(kind='value', lims=(-100, -75, -50))) -snr_stc.plot(figure=1, **kwargs) +kwargs = dict(initial_time=evoked.times[maxidx], hemi='split', + views=['lat', 'med'], subjects_dir=subjects_dir, size=(600, 600), + clim=dict(kind='value', lims=(-100, -70, -40)), + transparent=True, colormap='viridis') +snr_stc.plot(**kwargs) ############################################################################### -# EEG: +# EEG +# --- +# Next we do the same for EEG and plot the result on the cortex: evoked_eeg = evoked.copy().pick_types(eeg=True, meg=False) - -# Make inverse oparator: inv_op_eeg = make_inverse_operator(evoked_eeg.info, fwd, cov, fixed=True, verbose=True) - -# Calculate MNE: stc_eeg = apply_inverse(evoked_eeg, inv_op_eeg, lambda2, 'MNE', verbose=True) - -# Calculate SNR in source space: snr_stc_eeg = stc_eeg.estimate_snr(evoked_eeg.info, fwd, cov) - -# Plot an average SNR across source points over time: -ave_eeg = np.mean(snr_stc_eeg.data, axis=0) -fig, ax = plt.subplots() -ax.plot(evoked_eeg.times, ave_eeg) -ax.set(xlabel='Time (sec)', ylabel='SNR EEG') -fig.tight_layout() - -# Plot SNR on source space at the time point of maximum SNR (same time point as -# for MEG-EEG) -kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', - subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', - clim=dict(kind='value', lims=(-100, -75, -50))) -snr_stc_eeg.plot(figure=2, **kwargs) +snr_stc_eeg.plot(**kwargs) ############################################################################### -# MEG: +# MEG +# --- +# Finally we do this for MEG: evoked_meg = evoked.copy().pick_types(eeg=False, meg=True) - -# Make inverse oparator: inv_op_meg = make_inverse_operator(evoked_meg.info, fwd, cov, fixed=True, verbose=True) - -# Calculate MNE: stc_meg = apply_inverse(evoked_meg, inv_op_meg, lambda2, 'MNE', verbose=True) - -# Calculate SNR in source space: snr_stc_meg = stc_meg.estimate_snr(evoked_meg.info, fwd, cov) - -# Plot an average SNR across source points over time: -ave_meg = np.mean(snr_stc_meg.data, axis=0) -fig, ax = plt.subplots() -ax.plot(evoked_meg.times, ave_meg) -ax.set(xlabel='Time (sec)', ylabel='SNR MEG') -fig.tight_layout() - -# Plot SNR on source space at the time point of maximum SNR (same time point as -# for MEG-EEG) -kwargs = dict(initial_time=evoked.times[maxidx], hemi='both', - subjects_dir=subjects_dir, size=(600, 600), colormap='viridis', - clim=dict(kind='value', lims=(-100, -75, -50))) -snr_stc_meg.plot(figure=3, **kwargs) +snr_stc_meg.plot(**kwargs) + +############################################################################## +# References +# ---------- +# .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., +# Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). Mapping the +# Signal-To-Noise-Ratios of Cortical Sources in Magnetoencephalography +# and Electroencephalography. Human Brain Mapping, 30(4), 1077–1086. +# doi:10.1002/hbm.20571 diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 43025f0d0d3..f452f7e3072 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1505,7 +1505,6 @@ def extract_label_time_course(self, labels, src, mode='mean_flip', return label_tc - @verbose def estimate_snr(self, info, fwd, cov, verbose=None): r"""Compute time-varying SNR in the source space. @@ -1515,7 +1514,7 @@ def estimate_snr(self, info, fwd, cov, verbose=None): .. warning:: This function currently only works properly for fixed orientation. - + Parameters ---------- info : instance Info @@ -1545,10 +1544,10 @@ def estimate_snr(self, info, fwd, cov, verbose=None): forward model for a source with unit amplitude, :math:`a` is the source amplitude, :math:`N` is the number of sensors, and :math:`s_k^2` is the noise variance on sensor :math:`k`. - + References ---------- - .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, + .. [1] Goldenholz, D. M., Ahlfors, S. P., Hämäläinen, M. S., Sharon, D., Ishitobi, M., Vaina, L. M., & Stufflebeam, S. M. (2009). Mapping the Signal-To-Noise-Ratios of Cortical Sources in Magnetoencephalography and Electroencephalography. @@ -1556,28 +1555,25 @@ def estimate_snr(self, info, fwd, cov, verbose=None): """ from .forward import convert_forward_solution from .minimum_norm.inverse import _prepare_forward - + if (self.data >= 0).all(): - warn_('This STC appears to be from free orientation, currently SNR ' - 'function is valid only for fixed orientation') - - fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=True) - + warn_('This STC appears to be from free orientation, currently SNR' + ' function is valid only for fixed orientation') + + fwd = convert_forward_solution(fwd, surf_ori=True, force_fixed=False) + # G is gain matrix [ch x src], cov is noise covariance [ch x ch] - G, _, _, _, _, _, _, cov, _ = _prepare_forward(fwd, info, cov, - fixed=True, loose=0, rank=None, pca=False, - use_cps=True, exp=None, limit_depth_chs=False, - combine_xyz='fro', allow_fixed_depth=True, - limit=None) + G, _, _, _, _, _, _, cov, _ = _prepare_forward( + fwd, info, cov, fixed=True, loose=0, rank=None, pca=False, + use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', + allow_fixed_depth=False, limit=None) G = G['sol']['data'] - n_channels = cov['dim'] # number of sensors/channels - print('Number of channels in calculating SNR: %d' % n_channels) + n_channels = cov['dim'] # number of sensors/channels b_k2 = (G * G).T - s_k2 = np.diag(cov['data']) + s_k2 = np.diag(cov['data']) scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1, keepdims=True) - snr_stc = SourceEstimate(10 * np.log10((self.data * self.data) * scaling), - self.vertices, self.tmin, self.tstep, self.subject) - + snr_stc = self.copy() + snr_stc._data[:] = 10 * np.log10((self.data * self.data) * scaling) return snr_stc @@ -1714,7 +1710,6 @@ def center_of_mass(self, subject=None, hemi=None, restrict_vertices=False, return vertex, hemi, t - class _BaseVectorSourceEstimate(_BaseSourceEstimate): _data_ndim = 3 From 40c47134eb1b438c9c08f063d7542d16b024a1f8 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 4 Oct 2019 13:16:13 -0400 Subject: [PATCH 23/24] FIX: Rename --- .../{plot_source_space_SNR.py => plot_source_space_snr.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/inverse/{plot_source_space_SNR.py => plot_source_space_snr.py} (100%) diff --git a/examples/inverse/plot_source_space_SNR.py b/examples/inverse/plot_source_space_snr.py similarity index 100% rename from examples/inverse/plot_source_space_SNR.py rename to examples/inverse/plot_source_space_snr.py From 8c65e2eb8e4f738e109f37ef4ad65507966754c7 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 4 Oct 2019 14:00:10 -0400 Subject: [PATCH 24/24] FIX? --- examples/inverse/plot_source_space_snr.py | 2 -- mne/source_estimate.py | 3 +-- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/examples/inverse/plot_source_space_snr.py b/examples/inverse/plot_source_space_snr.py index 9813aca28ef..4191665d597 100644 --- a/examples/inverse/plot_source_space_snr.py +++ b/examples/inverse/plot_source_space_snr.py @@ -1,7 +1,5 @@ # -*- coding: utf-8 -*- """ -.. ex-mne-sourcespace-snr: - =============================== Computing source space SNR =============================== diff --git a/mne/source_estimate.py b/mne/source_estimate.py index f452f7e3072..6b71268aacf 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1568,7 +1568,7 @@ def estimate_snr(self, info, fwd, cov, verbose=None): use_cps=True, exp=None, limit_depth_chs=False, combine_xyz='fro', allow_fixed_depth=False, limit=None) G = G['sol']['data'] - n_channels = cov['dim'] # number of sensors/channels + n_channels = cov['dim'] # number of sensors/channels b_k2 = (G * G).T s_k2 = np.diag(cov['data']) scaling = (1 / n_channels) * np.sum(b_k2 / s_k2, axis=1, keepdims=True) @@ -1576,7 +1576,6 @@ def estimate_snr(self, info, fwd, cov, verbose=None): snr_stc._data[:] = 10 * np.log10((self.data * self.data) * scaling) return snr_stc - def get_peak(self, hemi=None, tmin=None, tmax=None, mode='abs', vert_as_index=False, time_as_index=False): """Get location and latency of peak amplitude.