diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 9f2d6afb7ff..262ea4e63c1 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -34,6 +34,8 @@ Changelog - Add option to :func:`mne.connectivity.spectral_connectivity` to compute corrected imaginary PLV by `Adonay Nunes`_ +- Add :func:`mne.SourceEstimate.estimate_snr` to estimate source-space SNR, by `Kaisu Lankinen`_ and `Padma Sundaram`_ + - Add option to specify the coordinate frame in :func:`mne.channels.read_custom_montage` by `Eric Larson`_ - Add reader for NIRx data in :func:`mne.io.read_raw_nirx` by `Robert Luke`_ diff --git a/doc/changes/names.inc b/doc/changes/names.inc index fa4e8a3ca3b..a79ed3cbccd 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -260,6 +260,10 @@ .. _Abram Hindle: http://softwareprocess.es +.. _Kaisu Lankinen: http://bishoplab.berkeley.edu/Kaisu.html + +.. _Padma Sundaram: https://www.nmr.mgh.harvard.edu/user/8071 + .. _Robert Luke: https://github.com/rob-luke .. _Mohammad Daneshzand: https://github.com/mdaneshzand diff --git a/mne/forward/forward.py b/mne/forward/forward.py index b8000eb85c9..18e15471252 100644 --- a/mne/forward/forward.py +++ b/mne/forward/forward.py @@ -43,7 +43,8 @@ write_trans) from ..utils import (_check_fname, get_subjects_dir, has_mne_c, warn, run_subprocess, check_fname, logger, verbose, fill_doc, - _validate_type, _check_compensation_grade, _check_option) + _validate_type, _check_compensation_grade, _check_option, + _check_stc_units) from ..label import Label from ..fixes import einsum @@ -1303,13 +1304,7 @@ def _apply_forward(fwd, stc, start=None, stop=None, on_missing='raise', 'Use pick_ori="normal" when computing the inverse to compute ' 'currents not current magnitudes.') - max_cur = np.max(np.abs(stc.data)) - if max_cur > 1e-7: # 100 nAm threshold for warning - warn('The maximum current magnitude is %0.1f nAm, which is very large.' - ' Are you trying to apply the forward model to noise-normalized ' - '(dSPM, sLORETA, or eLORETA) values? The result will only be ' - 'correct if currents (in units of Am) are used.' - % (1e9 * max_cur)) + _check_stc_units(stc) src_sel, stc_sel, _ = _stc_src_sel(fwd['src'], stc, on_missing=on_missing) gain = fwd['sol']['data'][:, src_sel] diff --git a/mne/source_estimate.py b/mne/source_estimate.py index d974d42cc13..1eb4c1deb4d 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -12,18 +12,21 @@ from scipy import linalg, sparse from scipy.sparse import coo_matrix, block_diag as sparse_block_diag +from .cov import Covariance +from .evoked import _get_peak from .filter import resample from .fixes import einsum -from .evoked import _get_peak from .surface import read_surface, _get_ico_surface, mesh_edges from .source_space import (_ensure_src, _get_morph_src_reordering, _ensure_src_subject, SourceSpaces) from .utils import (get_subjects_dir, _check_subject, logger, verbose, _time_mask, warn as warn_, copy_function_doc_to_method_doc, - fill_doc, _check_option, _validate_type, _check_src_normal) + fill_doc, _check_option, _validate_type, _check_src_normal, + _check_stc_units) from .viz import (plot_source_estimates, plot_vector_source_estimates, plot_volume_source_estimates) from .io.base import ToDataFrameMixin, TimeMixin +from .io.meas_info import Info from .externals.h5io import read_hdf5, write_hdf5 @@ -1553,9 +1556,12 @@ def estimate_snr(self, info, fwd, cov, verbose=None): Magnetoencephalography and Electroencephalography. Human Brain Mapping, 30(4), 1077–1086. doi:10.1002/hbm.20571 """ - from .forward import convert_forward_solution + from .forward import convert_forward_solution, Forward from .minimum_norm.inverse import _prepare_forward - + _validate_type(fwd, Forward, 'fwd') + _validate_type(info, Info, 'info') + _validate_type(cov, Covariance, 'cov') + _check_stc_units(self) if (self.data >= 0).all(): warn_('This STC appears to be from free orientation, currently SNR' ' function is valid only for fixed orientation') diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 39bc04b078d..737790e4171 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -12,8 +12,8 @@ read_evokeds, MixedSourceEstimate, find_events, Epochs, read_source_estimate, extract_label_time_course, spatio_temporal_tris_connectivity, - spatio_temporal_src_connectivity, - spatial_inter_hemi_connectivity, + spatio_temporal_src_connectivity, read_cov, + spatial_inter_hemi_connectivity, read_forward_solution, spatial_src_connectivity, spatial_tris_connectivity, SourceSpaces, VolVectorSourceEstimate) from mne.datasets import testing @@ -30,6 +30,13 @@ subjects_dir = op.join(data_path, 'subjects') fname_inv = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif') +fname_inv_fixed = op.join( + data_path, 'MEG', 'sample', + 'sample_audvis_trunc-meg-eeg-oct-4-meg-fixed-inv.fif') +fname_fwd = op.join( + data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg-eeg-oct-4-fwd.fif') +fname_cov = op.join( + data_path, 'MEG', 'sample', 'sample_audvis_trunc-cov.fif') fname_evoked = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-ave.fif') fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif') @@ -239,17 +246,25 @@ def _fake_vec_stc(n_time=10): 'foo') -def _real_vec_stc(): - inv = read_inverse_operator(fname_inv) +@testing.requires_testing_data +def test_stc_snr(): + """Test computing SNR from a STC.""" + inv = read_inverse_operator(fname_inv_fixed) + fwd = read_forward_solution(fname_fwd) + cov = read_cov(fname_cov) evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0].crop(0, 0.01) - return apply_inverse(evoked, inv, pick_ori='vector') - - -def _test_stc_integrety(stc): - """Test consistency of tmin, tstep, data.shape[-1] and times.""" - n_times = len(stc.times) - assert_equal(stc._data.shape[-1], n_times) - assert_array_equal(stc.times, stc.tmin + np.arange(n_times) * stc.tstep) + stc = apply_inverse(evoked, inv) + assert (stc.data < 0).any() + with pytest.warns(RuntimeWarning, match='nAm'): + stc.estimate_snr(evoked.info, fwd, cov) # dSPM + with pytest.warns(RuntimeWarning, match='free ori'): + abs(stc).estimate_snr(evoked.info, fwd, cov) + stc = apply_inverse(evoked, inv, method='MNE') + snr = stc.estimate_snr(evoked.info, fwd, cov) + assert_allclose(snr.times, evoked.times) + snr = snr.data + assert snr.max() < -10 + assert snr.min() > -120 def test_stc_attributes(): @@ -257,7 +272,9 @@ def test_stc_attributes(): stc = _fake_stc(n_time=10) vec_stc = _fake_vec_stc(n_time=10) - _test_stc_integrety(stc) + n_times = len(stc.times) + assert_equal(stc._data.shape[-1], n_times) + assert_array_equal(stc.times, stc.tmin + np.arange(n_times) * stc.tstep) assert_array_almost_equal( stc.times, [0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]) diff --git a/mne/utils/__init__.py b/mne/utils/__init__.py index 08aa68d6bbf..e0dc98088c9 100644 --- a/mne/utils/__init__.py +++ b/mne/utils/__init__.py @@ -14,7 +14,7 @@ _validate_type, _check_info_inv, _check_pylsl_installed, _check_channels_spatial_filter, _check_one_ch_type, _check_rank, _check_option, _check_depth, _check_combine, - _check_path_like, _check_src_normal) + _check_path_like, _check_src_normal, _check_stc_units) from .config import (set_config, get_config, get_config_path, set_cache_dir, set_memmap_min_size, get_subjects_dir, _get_stim_channel, sys_info, _get_extra_data_path, _get_root_dir, diff --git a/mne/utils/check.py b/mne/utils/check.py index 8d752fceb55..86f8d8164d4 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -537,3 +537,13 @@ def _check_src_normal(pick_ori, src): raise RuntimeError('Normal source orientation is supported only for ' 'surface or discrete SourceSpaces, got type ' '%s' % (src.kind,)) + + +def _check_stc_units(stc, threshold=1e-7): # 100 nAm threshold for warning + max_cur = np.max(np.abs(stc.data)) + if max_cur > threshold: + warn('The maximum current magnitude is %0.1f nAm, which is very large.' + ' Are you trying to apply the forward model to noise-normalized ' + '(dSPM, sLORETA, or eLORETA) values? The result will only be ' + 'correct if currents (in units of Am) are used.' + % (1e9 * max_cur))