MRG: Add convenience function to "prepare" emptyroom raw file for maxwell_filter

doc/changes/latest.inc

@@ -25,6 +25,8 @@ Enhancements
 ~~~~~~~~~~~~
 - Add support for ``overview_mode`` in :meth:`raw.plot() <mne.io.Raw.plot>` and related functions/methods (:gh:`10501` by `Eric Larson`_)
 
+- 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
 ~~~~
 - Fix bug in :func:`mne.io.read_raw_brainvision` when BrainVision data are acquired with the Brain Products "V-Amp" amplifier and disabled lowpass filter is marked with value ``0`` (:gh:`10517` by :newcontrib:`Alessandro Tonin`)

doc/preprocessing.rst

@@ -92,6 +92,7 @@ Projections:
    infomax
    equalize_bads
    maxwell_filter
+   maxwell_filter_prepare_emptyroom
    oversampled_temporal_projection
    peak_finder
    read_ica

mne/preprocessing/__init__.py

@@ -20,7 +20,7 @@
 from .infomax_ import infomax
 from .stim import fix_stim_artifact
 from .maxwell import (maxwell_filter, find_bad_channels_maxwell,
-                      compute_maxwell_basis)
+                      compute_maxwell_basis, maxwell_filter_prepare_emptyroom)
 from .realign import realign_raw
 from .xdawn import Xdawn
 from ._csd import compute_current_source_density

mne/preprocessing/maxwell.py

@@ -41,6 +41,88 @@
 # differences between algorithms
 
 
+@verbose
+def maxwell_filter_prepare_emptyroom(
+    raw_er, *, raw, bads='from_raw', verbose=None
+):
+    """Prepare an empty-room recording for Maxwell filtering.
+
+    Empty-room data by default lacks certain properties that are required to
+    ensure running :func:`~mne.preprocessing.maxwell_filter` will process the
+    empty-room recording the same way as the experimental data. This function
+    preconditions empty-room raw data instance accordingly so it can be used
+    for Maxwell filtering.
+
+    Specifically, this function will:
+
+    * Compile the list of bad channels according to the ``bads`` parameter.
+    * Inject the device-to-head transformation matrix from the experimental
+      recording into the empty-room recording.
+    * Set the same montage on the empty-room recording as the experimental
+      recording.
+
+    Parameters
+    ----------
+    raw_er : instance of Raw
+        The empty-room recording. It will not be modified.
+    raw : instance of Raw
+        The experimental recording, typically this will be the reference run
+        used for Maxwell filtering.
+    bads : 'from_raw' | 'union' | list of str
+        How to populate the list of bad channel names to be injected into
+        the empty-room recording. If ``'from_raw'`` (default) the list of bad
+        channels will be overwritten with that of ``raw``. If ``'union'``, will
+        use the union of bad channels in ``raw`` and ``raw_er``. Note that
+        this may lead to additional bad channels in the empty-room in
+        comparison to the experimental recording. If a list, you can explicitly
+        provide the list of channel names to be marked as bad.
+
+        .. note::
+           Non-MEG channels are silently dropped from the list of bads.
+    %(verbose)s
+
+    Returns
+    -------
+    raw_er_prepared : instance of Raw
+        A copy of the passed empty-room recording, ready for Maxwell filtering.
+
+    Notes
+    -----
+    .. versionadded:: 1.1
+    """
+    _validate_type(item=raw_er, types=BaseRaw, item_name='raw_er')
+    _validate_type(item=raw, types=BaseRaw, item_name='raw')
+    _validate_type(item=bads, types=(list, str), item_name='bads')
+    if isinstance(bads, str):
+        _check_option(
+            parameter='bads', value=bads, allowed_values=['from_raw', 'union']
+        )
+
+    raw_er_prepared = raw_er.copy()
+    del raw_er  # just to be sure, and to save memory
+
+    # handle bads; only keep MEG channels
+
+    if bads == 'from_raw':
+        bads = raw.info['bads']
+    elif bads == 'union':
+        bads = sorted(
+            set(raw.info['bads'] + raw_er_prepared.info['bads'])
+        )
+    bads = [ch_name for ch_name in bads
+            if ch_name.startswith('MEG')]
+    raw_er_prepared.info['bads'] = bads
+
+    # handle dev_head_t
+    raw_er_prepared.info['dev_head_t'] = raw.info['dev_head_t']
+
+    # handle montage
+    montage = raw.get_montage()
+    raw_er_prepared.set_montage(montage)
+
+    return raw_er_prepared
+
+
 # Changes to arguments here should also be made in find_bad_channels_maxwell
 @verbose
 def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,

mne/preprocessing/tests/test_maxwell.py

@@ -24,7 +24,8 @@
                     BaseRaw, read_raw_ctf)
 from mne.io.constants import FIFF
 from mne.preprocessing import (maxwell_filter, find_bad_channels_maxwell,
-                               annotate_amplitude, compute_maxwell_basis)
+                               annotate_amplitude, compute_maxwell_basis,
+                               maxwell_filter_prepare_emptyroom)
 from mne.preprocessing.maxwell import (
     _get_n_moments, _sss_basis_basic, _sh_complex_to_real,
     _sh_real_to_complex, _sh_negate, _bases_complex_to_real, _trans_sss_basis,
@@ -1443,3 +1444,46 @@ def test_compute_maxwell_basis(regularize, n):
     xform = S[:, :n_use_in] @ pS[:n_use_in]
     got = xform @ raw.pick_types(meg=True, exclude='bads').get_data()
     assert_allclose(got, want)
+
+
+@testing.requires_testing_data
+@pytest.mark.parametrize(
+    'bads', ['from_raw', 'union', ['MEG 0113', 'MEG 2313']]
+)
+def test_prepare_emptyroom(bads):
+    """Test prepare_emptyroom."""
+    raw = read_raw_fif(sample_fname)
+    raw_er = raw.copy().pick_types(meg=True)
+    raw_er.info['dev_head_t'] = None
+    raw_er.set_montage(None)
+
+    if bads == 'from_raw':
+        raw_bads_orig = ['MEG 0113', 'MEG 2313']
+        raw_er_bads_orig = []
+    elif bads == 'union':
+        raw_bads_orig = ['MEG 0113']
+        raw_er_bads_orig = ['MEG 2313']
+    else:
+        raw_bads_orig = bads
+        raw_er_bads_orig = []
+
+    raw.info['bads'] = raw_bads_orig
+    raw_er.info['bads'] = raw_er_bads_orig
+
+    raw_er_prepared = maxwell_filter_prepare_emptyroom(
+        raw_er=raw_er,
+        raw=raw,
+        bads=bads
+    )
+    assert raw_er_prepared.info['bads'] == ['MEG 0113', 'MEG 2313']
+    assert raw_er_prepared.info['dev_head_t'] == raw.info['dev_head_t']
+
+    montage_expected = raw.copy().pick_types(meg=True).get_montage()
+    # XXX assert raw_er.get_montage == montage_expected fails – why?
+    assert raw_er_prepared.get_montage().dig == montage_expected.dig
+
+    # Ensure the originals were not modified
+    assert raw.info['bads'] == raw_bads_orig
+    assert raw_er.info['bads'] == raw_er_bads_orig
+    assert raw_er.info['dev_head_t'] is None
+    assert raw_er.get_montage() is None

setup.cfg

@@ -15,7 +15,7 @@ release = egg_info -RDb ''
 doc_files = doc
 
 [flake8]
-exclude = __init__.py,constants.py,fixes.py,resources.py,*doc/auto_examples*,*doc/_build*
+exclude = __init__.py,constants.py,fixes.py,resources.py,*doc/auto_examples*,*doc/_build*,logo/*
 ignore = W503,W504,I100,I101,I201,N806,E201,E202,E221,E222,E241
 # We add A for the array-spacing plugin, and ignore the E ones it covers above
 select = A,E,F,W,C