Merge remote-tracking branch 'upstream/main' into hoechenberger/issue10533

* upstream/main:
  MAINT: Extra test for coreg (mne-tools#10549)
  BUG: Fix annot meas_date / crop (mne-tools#10491)
  Update latest.inc

Author: larsoner

doc/changes/latest.inc

@@ -27,6 +27,12 @@ Enhancements
 
 - The ``pick_channels`` method gained a ``verbose`` parameter, allowing e.g. to suppress messages about removed projectors (:gh:`10544` by `Richard Höchenberger`_)
 
+- The :func:`mne.make_forward_dipole` function can now take a list of dipoles to make a multi-dipole forward models (:gh:`10464` by `Marijn van Vliet`_)
+
+- Add example of Xfit-style ECD modeling using multiple dipoles (:gh:`10464` by `Marijn van Vliet`_)
+
+- It is now possible to compute inverse solutions with restricted source orientations using discrete forward models (:gh:`10464` by `Marijn van Vliet`_)
+
 - The new function :func:`mne.preprocessing.maxwell_filter_prepare_emptyroom` simplifies the preconditioning of an empty-room recording for our Maxwell filtering operations (:gh:`10533` by `Richard Höchenberger`_)
 
 Bugs
@@ -39,6 +45,8 @@ Bugs
 
 - Fix bug where ``theme`` was not handled properly in :meth:`mne.io.Raw.plot` (:gh:`10487`, :gh:`10500` by `Mathieu Scheltienne`_ and `Eric Larson`_)
 
+- Fix bug in :meth:`raw.crop(start, stop) <mne.io.Raw.crop>` that would cause annotations to be erroneously shifted when ``start != 0`` and no measurement date was set. (:gh:`10491` by `Eric Larson`_)
+
 - Rendering issues with recent MESA releases can be avoided by setting the new environment variable``MNE_3D_OPTION_MULTI_SAMPLES=1`` or using :func:`mne.viz.set_3d_options` (:gh:`10513` by `Eric Larson`_)
 
 - Fix behavior for the ``pyvista`` 3D renderer's ``quiver3D`` function so that default arguments plot a glyph in ``arrow`` mode (:gh:`10493` by `Alex Rockhill`_)

mne/annotations.py

@@ -111,7 +111,7 @@ class Annotations(object):
         starting time of annotation acquisition. If None (default),
         starting time is determined from beginning of raw data acquisition.
         In general, ``raw.info['meas_date']`` (or None) can be used for syncing
-        the annotations with raw data if their acquisiton is started at the
+        the annotations with raw data if their acquisition is started at the
         same time. If it is a string, it should conform to the ISO8601 format.
         More precisely to this '%%Y-%%m-%%d %%H:%%M:%%S.%%f' particular case of
         the ISO8601 format where the delimiter between date and time is ' '.
@@ -231,6 +231,14 @@ class Annotations(object):
              e                        +------+
          orig_time                 onset[0]'
 
+    .. warning::
+       This means that when ``raw.info['meas_date'] is None``, doing
+       ``raw.set_annotations(raw.annotations)`` will not alter ``raw`` if and
+       only if ``raw.first_samp == 0``. When it's non-zero,
+       ``raw.set_annotations`` will assume that the "new" annotations refer to
+       the original data (with ``first_samp==0``), and will be re-referenced to
+       the new time offset!
+
     **Specific annotation**
 
     ``BAD_ACQ_SKIP`` annotation leads to specific reading/writing file
@@ -1187,7 +1195,7 @@ def _read_brainstorm_annotations(fname, orig_time=None):
         starting time of annotation acquisition. If None (default),
         starting time is determined from beginning of raw data acquisition.
         In general, ``raw.info['meas_date']`` (or None) can be used for syncing
-        the annotations with raw data if their acquisiton is started at the
+        the annotations with raw data if their acquisition is started at the
         same time.
 
     Returns

mne/io/base.py

@@ -1353,9 +1353,15 @@ def crop(self, tmin=0.0, tmax=None, include_tmax=True, *, verbose=None):
             self._data = self._data[:, smin:smax + 1].copy()
 
         annotations = self.annotations
-        if annotations.orig_time is None:
-            annotations.onset -= tmin
         # now call setter to filter out annotations outside of interval
+        if annotations.orig_time is None:
+            assert self.info['meas_date'] is None
+            # When self.info['meas_date'] is None (which is guaranteed if
+            # self.annotations.orig_time is None), when we do the
+            # self.set_annotations, it's assumed that the annotations onset
+            # are relative to first_time, so we have to subtract it, then
+            # set_annotations will put it back.
+            annotations.onset -= self.first_time
         self.set_annotations(annotations, False)
 
         return self

mne/io/ctf/info.py

@@ -98,7 +98,7 @@ def _convert_time(date_str, time_str):
         raise RuntimeError('Illegal time: %s' % time_str)
     # MNE-C uses mktime which uses local time, but here we instead decouple
     # conversion location from the process, and instead assume that the
-    # acquisiton was in GMT. This will be wrong for most sites, but at least
+    # acquisition was in GMT. This will be wrong for most sites, but at least
     # the value we obtain here won't depend on the geographical location
     # that the file was converted.
     res = timegm((date.tm_year, date.tm_mon, date.tm_mday,

mne/io/meas_info.py

@@ -723,7 +723,7 @@ class Info(dict, MontageMixin, ContainsMixin):
         smartshield : dict
             MaxShield information. This dictionary is (always?) empty,
             but its presence implies that MaxShield was used during
-            acquisiton.
+            acquisition.
 
     * ``subject_info`` dict:
 

mne/tests/test_annotations.py

@@ -1440,22 +1440,34 @@ def test_annotation_duration_setting():
 
 
 @pytest.mark.parametrize('meas_date', (None, 1))
-@pytest.mark.parametrize('first_samp', (0, 100))
-def test_annot_noop(meas_date, first_samp):
+@pytest.mark.parametrize('set_meas_date', ('before', 'after'))
+@pytest.mark.parametrize('first_samp', (0, 100, 3000))
+def test_annot_noop(meas_date, first_samp, set_meas_date):
     """Show some unintuitive behavior of annotations."""
     sfreq = 1000.
-    raw = RawArray(np.zeros((1, 2000)), create_info(1, sfreq, 'eeg'),
-                   first_samp=first_samp)
-    annot = Annotations(0.5, 0.1, 'bad')
-    raw.set_annotations(annot)
-    raw.set_meas_date(meas_date)
+    info = create_info(1, sfreq, 'eeg')
+    onset = 0.5
+    annot_kwargs = dict()
+    if set_meas_date == 'before':
+        with info._unlock():
+            info['meas_date'] = _handle_meas_date(meas_date)
+        if meas_date is not None:
+            onset += first_samp / sfreq
+        annot_kwargs['orig_time'] = meas_date
+    raw = RawArray(np.zeros((1, 2000)), info, first_samp=first_samp)
+    annot = Annotations(onset, 0.1, 'bad', **annot_kwargs)
+    raw.set_annotations(annot, verbose='debug')
+    if set_meas_date == 'after':
+        raw.set_meas_date(meas_date)
     first_annot = raw.annotations
-    raw.set_annotations(first_annot)  # should be a no-op...
+    if meas_date is None:
+        first_annot.onset -= raw.first_time
+    raw.set_annotations(first_annot, verbose='debug')  # should be a no-op...
     second_annot = raw.annotations
     want = first_annot.onset[0]
     # it has been shifted when meas_date is None!
     if meas_date is None:
-        want = want + first_samp / sfreq
+        want = want + raw.first_time
     assert_allclose(second_annot.onset[0], want)
 
 
@@ -1532,9 +1544,8 @@ def _create_raw(eeg, sfreq, onset, description, meas_date, first_samp,
     end_idx = np.where(raw.annotations.description == 'off')[0]
     tmins = raw.annotations.onset[start_idx]
     tmaxs = raw.annotations.onset[end_idx]
-    if meas_date_1 is not None:  # need to reference to what will be 0 time
-        tmins -= raw.first_time
-        tmaxs -= raw.first_time
+    tmins -= raw.first_time
+    tmaxs -= raw.first_time
     assert len(tmins) == len(tmaxs) == 2
     assert raw.info['meas_date'] == meas_date_1
     _assert_annotations_equal(raw.annotations, want_annot)
@@ -1571,3 +1582,41 @@ def _create_raw(eeg, sfreq, onset, description, meas_date, first_samp,
         assert sess.first_samp == want_first_samp
         assert sess.annotations.orig_time == meas_date_1
         assert list(sess.annotations.description) == descs
+
+
+@pytest.mark.parametrize('first_samp', (0, 10000))
+@pytest.mark.parametrize('meas_date', (None, 24 * 60 * 60))
+def test_annot_meas_date_first_samp_crop(meas_date, first_samp):
+    """Test yet another meas_date / first_samp issue."""
+    sfreq = 1000.
+    info = mne.create_info(1, sfreq, 'eeg')
+    raw = mne.io.RawArray(
+        np.random.RandomState(0).randn(1, 3000), info, first_samp=first_samp)
+    raw.set_meas_date(meas_date)
+    onset = np.array([0, 1, 2], float)
+    if meas_date is not None:
+        onset += first_samp / sfreq
+    annot = mne.Annotations(
+        onset=onset,
+        duration=[0.1, 0.2, 0.3],
+        description=["a", "b", "c"],
+        orig_time=raw.info['meas_date'])
+    assert len(annot) == 3
+    raw.set_annotations(annot)
+    assert len(raw.annotations) == 3
+    raw_crop = raw.copy().crop(0, 1.5, verbose='debug')
+    assert len(raw_crop.annotations) == 2
+    assert_array_equal(raw_crop.annotations.description, annot.description[:2])
+    assert_array_equal(raw_crop.annotations.duration, annot.duration[:2])
+    # these two should be the equivalent
+    raw_crop = raw.copy().crop(2, 2.5, verbose='debug')
+    raw_crop_2 = raw.copy().crop(1, None).crop(1, 1.5)
+    assert_allclose(raw_crop.get_data(), raw_crop_2.get_data())
+    assert raw_crop.first_samp == raw_crop_2.first_samp
+    want_onset = onset[2:]
+    if meas_date is None:
+        want_onset = want_onset + raw.first_time
+    for this_raw in (raw_crop, raw_crop_2):
+        assert len(this_raw.annotations) == 1
+        assert_allclose(this_raw.annotations.onset, want_onset)
+        assert_allclose(this_raw.annotations.duration, annot.duration[2:])

mne/tests/test_coreg.py

@@ -13,7 +13,7 @@
 from mne.datasets import testing
 from mne.transforms import (Transform, apply_trans, rotation, translation,
                             scaling, read_trans, _angle_between_quats,
-                            rot_to_quat)
+                            rot_to_quat, invert_transform)
 from mne.coreg import (fit_matched_points, create_default_subject, scale_mri,
                        _is_mri_subject, scale_labels, scale_source_space,
                        coregister_fiducials, get_mni_fiducials, Coregistration)
@@ -235,6 +235,17 @@ def test_scale_mri_xfm(tmp_path, few_surfaces):
         mni = mne.vertex_to_mni(vertices, hemis, subject_to,
                                 subjects_dir=tempdir)
         assert_allclose(mni, mni_from, atol=1e-3)  # 0.001 mm
+        # Check head_to_mni (the `trans` here does not really matter)
+        trans = rotation(0.001, 0.002, 0.003) @ translation(0.01, 0.02, 0.03)
+        trans = Transform('head', 'mri', trans)
+        pos_head_from = np.random.RandomState(0).randn(4, 3)
+        pos_mni_from = mne.head_to_mni(
+            pos_head_from, subject_from, trans, tempdir)
+        pos_mri_from = apply_trans(trans, pos_head_from)
+        pos_mri = pos_mri_from * scale
+        pos_head = apply_trans(invert_transform(trans), pos_mri)
+        pos_mni = mne.head_to_mni(pos_head, subject_to, trans, tempdir)
+        assert_allclose(pos_mni, pos_mni_from, atol=1e-3)
 
 
 def test_fit_matched_points():