diff --git a/doc/changes/dev/14180.bugfix.rst b/doc/changes/dev/14180.bugfix.rst new file mode 100644 index 00000000000..2d9cd720e5f --- /dev/null +++ b/doc/changes/dev/14180.bugfix.rst @@ -0,0 +1 @@ +Fix a bug that caused :meth:`mne.io.Raw.pick` and :meth:`mne.io.Raw.drop_channels` to fail on the result of :func:`mne.preprocessing.maxwell_filter` when movement compensation was enabled with ``head_pos``, by `Christian Brodbeck`_. diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index b161e07d674..c5d9245d527 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -312,6 +312,11 @@ def maxwell_filter( ----- .. versionadded:: 0.11 + When ``head_pos`` is provided, the returned object contains cHPI result + channels that have no counterpart in the source file. It therefore cannot + be concatenated using ``raw_sss.append(..., preload=False)``. Use + ``preload=True`` or a memory-mapped filename instead. + Some of this code was adapted and relicensed (with BSD form) with permission from Jussi Nurminen. These algorithms are based on work from :footcite:`TauluKajola2005` and :footcite:`TauluSimola2006`. @@ -1258,6 +1263,20 @@ def _copy_preload_add_channels(raw, add_channels, copy, info): raw.info["chs"].extend(chpi_chs) raw.info._update_redundant() raw.info._check_consistency() + # The remaining per-channel attributes must grow along with info, otherwise + # any later channel operation (e.g., raw_sss.drop_channels) indexes them out + # of bounds + raw._cals = np.concatenate([raw._cals, raw.info._cals[off:]]) + # The added channels have no counterpart in the source file, but the data are + # preloaded, so _read_picks (and _raw_extras) will never be used to read from + # disk again -- use indices that are likely to break loudly if they ever are + extra_idx = [2147483647] * len(chpi_chs) # 2 ** 31 - 1 + raw._read_picks = [np.concatenate([r, extra_idx]) for r in raw._read_picks] + assert raw._comp is None # preloading the data above unsets it + if raw._projector is not None: # identity for the added channels + projector = np.eye(raw.info["nchan"]) + projector[:off, :off] = raw._projector + raw._projector = projector assert raw._data.shape == (raw.info["nchan"], len(raw.times)) # Return the pos picks pos_picks = np.arange(len(raw.ch_names) - len(chpi_chs), len(raw.ch_names)) diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py index 989ab780e30..55718d901a9 100644 --- a/mne/preprocessing/tests/test_maxwell.py +++ b/mne/preprocessing/tests/test_maxwell.py @@ -256,6 +256,16 @@ def test_movement_compensation_basic(tmp_path): assert_meg_snr( raw_sss, read_crop(sss_movecomp_fname, lims), 4.6, 12.4, chpi_med_tol=58 ) + # appending the head position channels must leave the instance consistent, i.e., + # picking and dropping channels must still work (and agree with info) + want = raw_sss.ch_names[:2] + ["CHPI002"] + raw_sss_pick = raw_sss.copy().pick(want) + assert raw_sss_pick.ch_names == want + assert len(raw_sss_pick._cals) == len(raw_sss_pick._read_picks[0]) == len(want) + assert_allclose(raw_sss_pick.get_data(), raw_sss.get_data(want)) + raw_sss_drop = raw_sss.copy().drop_channels(["CHPI000"]) + assert raw_sss_drop.ch_names == [c for c in raw_sss.ch_names if c != "CHPI000"] + assert_allclose(raw_sss_drop.get_data(), raw_sss.get_data(raw_sss_drop.ch_names)) # IO temp_fname = tmp_path / "test_raw_sss.fif" raw_sss.save(temp_fname)