MRG, FIX: Add tests for estimate_snr

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`_

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

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]

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')

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,25 +246,35 @@ 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():
     """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])
 

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,

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))