FIX: Add support for surfaces, add test

mne/label.py

@@ -2158,7 +2158,9 @@ def labels_to_stc(labels, values, tmin=0, tstep=1, subject=None, src=None,
     subject : str | None
         The subject for which to create the STC.
     %(eltc_src)s
-        Only used if the values come from a volumetric source space.
+        Can be omitted if using a surface source space, in which case
+        the label vertices will determine the output STC vertices.
+        Required if using a volumetric source space.
 
         .. versionadded:: 0.22
     %(verbose)s
@@ -2185,25 +2187,27 @@ def labels_to_stc(labels, values, tmin=0, tstep=1, subject=None, src=None,
         raise ValueError('values must have 1 or 2 dimensions, got %s'
                          % (values.ndim,))
     _validate_type(src, (SourceSpaces, None))
-    kind = getattr(src, 'kind', 'surface')  # src can be None
-    _check_option('source space kind', kind, ('surface', 'volume'))
-    if kind == 'surface':
+    if src is None:
         data, vertices = _labels_to_stc_surf(
             labels, values, tmin, tstep, subject)
         klass = SourceEstimate
     else:
-        assert kind == 'volume'
+        kind = src.kind
+        _check_option('source space kind', kind, ('surface', 'volume'))
+        if kind == 'volume':
+            klass = VolSourceEstimate
+        else:
+            klass = SourceEstimate
         # Easiest way is to get a dot-able operator and use it
         vertices = [s['vertno'].copy() for s in src]
-        stc = VolSourceEstimate(
+        stc = klass(
             np.eye(sum(len(v) for v in vertices)), vertices, 0, 1, subject)
         label_op = extract_label_time_course(
-            stc, labels, src, allow_empty=True)
+            stc, labels, src=src, mode='mean', allow_empty=True)
         _check_values_labels(values, label_op.shape[0])
         rev_op = np.zeros(label_op.shape[::-1])
         rev_op[np.arange(label_op.shape[1]), np.argmax(label_op, axis=0)] = 1.
         data = rev_op @ values
-        klass = VolSourceEstimate
     return klass(data, vertices, tmin, tstep, subject, verbose)
 
 

mne/morph.py

@@ -220,7 +220,8 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage',
             src_data['to_vox_map'] = (src_to[-1]['shape'], src_ras_t)
             vertices_to_vol = [s['vertno'] for s in src_to[surf_offset:]]
             zooms_src_to = np.diag(src_to[-1]['src_mri_t']['trans'])[:3] * 1000
-            assert (zooms_src_to[0] == zooms_src_to).all()
+            if (zooms_src_to[0] != zooms_src_to).any():
+                raise RuntimeError('Non-uniform zooms not supported')
             zooms_src_to = tuple(zooms_src_to)
 
         # pre-compute non-linear morph

mne/tests/test_source_estimate.py

@@ -16,6 +16,7 @@
 from scipy.optimize import fmin_cobyla
 from scipy.spatial.distance import cdist
 
+import mne
 from mne import (stats, SourceEstimate, VectorSourceEstimate,
                  VolSourceEstimate, Label, read_source_spaces,
                  read_evokeds, MixedSourceEstimate, find_events, Epochs,
@@ -27,14 +28,16 @@
                  SourceSpaces, VolVectorSourceEstimate, read_trans, pick_types,
                  MixedVectorSourceEstimate, setup_volume_source_space,
                  convert_forward_solution, pick_types_forward,
-                 compute_source_morph, labels_to_stc)
+                 compute_source_morph, labels_to_stc, scale_mri,
+                 read_freesurfer_lut, write_source_spaces)
 from mne.datasets import testing
 from mne.externals.h5io import write_hdf5
 from mne.fixes import fft, _get_img_fdata, nullcontext
 from mne.io import read_info
 from mne.io.constants import FIFF
 from mne.source_estimate import grade_to_tris, _get_vol_mask
 from mne.source_space import _get_src_nn
+from mne.surface import _make_morph_map_hemi
 from mne.transforms import apply_trans, invert_transform, transform_surface_to
 from mne.minimum_norm import (read_inverse_operator, apply_inverse,
                               apply_inverse_epochs, make_inverse_operator)
@@ -1579,3 +1582,153 @@ def test_stc_near_sensors(tmpdir):
     assert isinstance(stc_vol, VolSourceEstimate)
     log = log.getvalue()
     assert '4157 volume vertices' in log
+
+
+def _make_morph_map_hemi_same(subject_from, subject_to, subjects_dir,
+                              reg_from, reg_to):
+    return _make_morph_map_hemi(subject_from, subject_from, subjects_dir,
+                                reg_from, reg_from)
+
+
+@pytest.mark.parametrize('kind', ('volume', 'surface'))
+@pytest.mark.parametrize('scale', (1., 0.75, (1.0, 0.8, 1.2)))
+def test_scale_morph_labels(kind, scale, monkeypatch, tmpdir):
+    """Test label extraction, morphing, and MRI scaling relationships."""
+    tempdir = str(tmpdir)
+    subject_from = 'sample'
+    subject_to = 'small'
+    testing_dir = op.join(subjects_dir, subject_from)
+    from_dir = op.join(tempdir, subject_from)
+    for root in ('mri', 'surf', 'label', 'bem'):
+        os.makedirs(op.join(from_dir, root), exist_ok=True)
+    for hemi in ('lh', 'rh'):
+        for root, fname in (('surf', 'sphere'), ('surf', 'white'),
+                            ('surf', 'sphere.reg'),
+                            ('label', 'aparc.annot')):
+            use_fname = op.join(root, f'{hemi}.{fname}')
+            copyfile(op.join(testing_dir, use_fname),
+                     op.join(from_dir, use_fname))
+    for root, fname in (('mri', 'aseg.mgz'), ('mri', 'brain.mgz'),
+                        ('bem', f'{subject_from}-oct-4-src.fif')):
+        use_fname = op.join(root, fname)
+        copyfile(op.join(testing_dir, use_fname),
+                 op.join(from_dir, use_fname))
+    del testing_dir
+    if kind == 'surface':
+        src_from = read_source_spaces(fname_src_3)
+        assert src_from[0]['dist'] is None
+        klass = SourceEstimate
+        assert len(src_from) == 2
+        assert src_from[0]['nuse'] == src_from[1]['nuse'] == 258
+        labels_from = read_labels_from_annot(
+            subject_from, subjects_dir=tempdir)
+        n_labels = len(labels_from)
+    else:
+        assert kind == 'volume'
+        volume_labels = list(read_freesurfer_lut()[0])
+        src_from = setup_volume_source_space(
+            subject_from, 20., 'aseg.mgz',
+            volume_label=volume_labels, single_volume=True,
+            subjects_dir=tempdir)
+        labels_from = op.join(
+            tempdir, subject_from, 'mri', 'aseg.mgz')
+        n_labels = 46
+        assert op.isfile(labels_from)
+        klass = VolSourceEstimate
+        assert len(src_from) == 1
+        assert src_from[0]['nuse'] == 389
+        write_source_spaces(
+            op.join(from_dir, 'bem', 'sample-vol20-src.fif'), src_from)
+    scale_mri(subject_from, subject_to, scale, subjects_dir=tempdir,
+              annot=True, skip_fiducials=True, verbose=True,
+              overwrite=True)  # XXX
+    if kind == 'surface':
+        src_to = read_source_spaces(
+            op.join(tempdir, subject_to, 'bem',
+                    f'{subject_to}-oct-4-src.fif'))
+        labels_to = read_labels_from_annot(
+            subject_to, subjects_dir=tempdir)
+        # Save time since we know these subjects are identical
+        monkeypatch.setattr(mne.surface, '_make_morph_map_hemi',
+                            _make_morph_map_hemi_same)
+    else:
+        src_to = read_source_spaces(
+            op.join(tempdir, subject_to, 'bem',
+                    f'{subject_to}-vol20-src.fif'))
+        labels_to = op.join(
+            tempdir, subject_to, 'mri', 'aseg.mgz')
+    # 1. Label->STC->Label for the given subject should be identity
+    #    (for surfaces at least; for volumes it's not as clean as this
+    #     due to interpolation)
+    n_times = 50
+    rng = np.random.RandomState(0)
+    label_tc = rng.randn(n_labels, n_times)
+    # check that a random permutation of our labels yields a terrible
+    # correlation
+    corr = np.corrcoef(label_tc.ravel(),
+                       rng.permutation(label_tc).ravel())[0, 1]
+    assert -0.06 < corr < 0.06
+    # project label activations to full source space
+    stc = labels_to_stc(labels_from, label_tc, src=src_from)
+    assert isinstance(stc, klass)
+    label_tc_from = extract_label_time_course(
+        stc, labels_from, src_from, mode='mean')
+    if kind == 'surface':
+        assert_allclose(label_tc, label_tc_from, rtol=1e-12, atol=1e-12)
+    else:
+        corr = np.corrcoef(label_tc.ravel(), label_tc_from.ravel())[0, 1]
+        assert 0.66 < corr < 0.69
+
+    #
+    # 2. Changing STC subject to the surrogate and then extracting
+    #
+    stc.subject = subject_to
+    label_tc_to = extract_label_time_course(
+        stc, labels_to, src_to, mode='mean')
+    assert_allclose(label_tc_from, label_tc_to, rtol=1e-12, atol=1e-12)
+    stc.subject = subject_from
+
+    #
+    # 3. Morphing STC to new subject then extracting
+    #
+    if isinstance(scale, tuple) and kind == 'volume':
+        ctx = pytest.raises(RuntimeError, match='Non-uniform zooms')
+        test_morph = False
+    elif kind == 'surface':
+        ctx = pytest.warns(RuntimeWarning, match='not included')
+        test_morph = True
+    else:
+        ctx = nullcontext()
+        test_morph = True
+    with ctx:  # vertices not included
+        morph = compute_source_morph(
+            src_from, subject_to=subject_to, src_to=src_to,
+            subjects_dir=tempdir, niter_sdr=(), smooth=1, verbose=True)
+    if test_morph:
+        stc_to = morph.apply(stc)
+        label_tc_to_morph = extract_label_time_course(
+            stc_to, labels_to, src_to, mode='mean')
+        if kind == 'volume':
+            corr = np.corrcoef(
+                label_tc.ravel(), label_tc_to_morph.ravel())[0, 1]
+            if scale == 1.:
+                assert 0.42 < corr <= 0.47, corr
+            else:
+                # XXX even just a pure scaling does not work very well, bug...
+                assert -0.1 < corr < 0.1
+        else:
+            assert_allclose(
+                label_tc, label_tc_to_morph, atol=1e-12, rtol=1e-12)
+
+    #
+    # 4. The same round trip from (1) but in the warped space
+    #
+    stc = labels_to_stc(labels_to, label_tc, src=src_to)
+    assert isinstance(stc, klass)
+    label_tc_to = extract_label_time_course(
+        stc, labels_to, src_to, mode='mean')
+    if kind == 'surface':
+        assert_allclose(label_tc, label_tc_to, rtol=1e-12, atol=1e-12)
+    else:
+        corr = np.corrcoef(label_tc.ravel(), label_tc_to.ravel())[0, 1]
+        assert 0.66 < corr < 0.69