BUG: Fix bug with recon trans

.mailmap

@@ -220,6 +220,7 @@ Mikołaj Magnuski <[email protected]> Mikolaj Magnuski <[email protected]
 Mikołaj Magnuski <[email protected]> mmagnuski <[email protected]>
 Mohamed Sherif <[email protected]> mohdsherif <[email protected]>
 Mohammad Daneshzand <[email protected]> mdaneshzand <[email protected]>
+Motofumi Fushimi <[email protected]> motofumi-fushimi <[email protected]>
 Natalie Klein <[email protected]> natalieklein <[email protected]>
 Nathalie Gayraud <[email protected]> Nathalie <[email protected]>
 Nathalie Gayraud <[email protected]> Nathalie <[email protected]>

doc/changes/devel/12348.bugfix.rst

@@ -0,0 +1 @@
+Fix bug in :func:`mne.preprocessing.maxwell_filter` where calibration was incorrectly applied during virtual sensor reconstruction, by `Eric Larson`_ and :newcontrib:`Motofumi Fushimi`.

doc/changes/names.inc

@@ -386,6 +386,8 @@
 
 .. _Moritz Gerster: https://github.com/moritz-gerster
 
+.. _Motofumi Fushimi: https://github.com/motofumi-fushimi/motofumi-fushimi.github.io
+
 .. _Natalie Klein: https://github.com/natalieklein
 
 .. _Nathalie Gayraud: https://github.com/ngayraud

mne/preprocessing/maxwell.py

@@ -519,8 +519,12 @@ def _prep_maxwell_filter(
     #
     sss_cal = dict()
     if calibration is not None:
+        # Modifies info in place, so make a copy for recon later
+        info_recon = info.copy()
         calibration, sss_cal = _update_sensor_geometry(info, calibration, ignore_ref)
         mag_or_fine.fill(True)  # all channels now have some mag-type data
+    else:
+        info_recon = info
 
     # Determine/check the origin of the expansion
     origin = _check_origin(origin, info, coord_frame, disp=True)
@@ -553,7 +557,8 @@ def _prep_maxwell_filter(
     #
     exp = dict(origin=origin_head, int_order=int_order, ext_order=0)
     all_coils = _prep_mf_coils(info, ignore_ref)
-    S_recon = _trans_sss_basis(exp, all_coils, recon_trans, coil_scale)
+    all_coils_recon = _prep_mf_coils(info_recon, ignore_ref)
+    S_recon = _trans_sss_basis(exp, all_coils_recon, recon_trans, coil_scale)
     exp["ext_order"] = ext_order
     exp["extended_proj"] = extended_proj
     del extended_proj

mne/preprocessing/tests/test_maxwell.py

@@ -730,7 +730,8 @@ def test_spatiotemporal_only():
     raw_tsss = maxwell_filter(raw, st_duration=tmax, st_correlation=1.0, st_only=True)
     assert_allclose(raw[:][0], raw_tsss[:][0])
     # degenerate
-    pytest.raises(ValueError, maxwell_filter, raw, st_only=True)  # no ST
+    with pytest.raises(ValueError, match="must not be None if st_only"):
+        maxwell_filter(raw, st_only=True)
     # two-step process equivalent to single-step process
     raw_tsss = maxwell_filter(raw, st_duration=tmax, st_only=True)
     raw_tsss = maxwell_filter(raw_tsss)
@@ -771,7 +772,7 @@ def test_fine_calibration():
     log = log.getvalue()
     assert "Using fine calibration" in log
     assert fine_cal_fname.stem in log
-    assert_meg_snr(raw_sss, sss_fine_cal, 82, 611)
+    assert_meg_snr(raw_sss, sss_fine_cal, 1.3, 180)  # similar to MaxFilter
     py_cal = raw_sss.info["proc_history"][0]["max_info"]["sss_cal"]
     assert py_cal is not None
     assert len(py_cal) > 0
@@ -812,15 +813,11 @@ def test_fine_calibration():
         regularize=None,
         bad_condition="ignore",
     )
-    assert_meg_snr(raw_sss_3D, sss_fine_cal, 1.0, 6.0)
+    assert_meg_snr(raw_sss_3D, sss_fine_cal, 0.9, 6.0)
+    assert_meg_snr(raw_sss_3D, raw_sss, 1.1, 6.0)  # slightly better than 1D
     raw_ctf = read_crop(fname_ctf_raw).apply_gradient_compensation(0)
-    pytest.raises(
-        RuntimeError,
-        maxwell_filter,
-        raw_ctf,
-        origin=(0.0, 0.0, 0.04),
-        calibration=fine_cal_fname,
-    )
+    with pytest.raises(RuntimeError, match="Not all MEG channels"):
+        maxwell_filter(raw_ctf, origin=(0.0, 0.0, 0.04), calibration=fine_cal_fname)
 
 
 @pytest.mark.slowtest
@@ -884,7 +881,8 @@ def test_cross_talk(tmp_path):
     assert len(py_ctc) > 0
     with pytest.raises(TypeError, match="path-like"):
         maxwell_filter(raw, cross_talk=raw)
-    pytest.raises(ValueError, maxwell_filter, raw, cross_talk=raw_fname)
+    with pytest.raises(ValueError, match="Invalid cross-talk FIF"):
+        maxwell_filter(raw, cross_talk=raw_fname)
     mf_ctc = sss_ctc.info["proc_history"][0]["max_info"]["sss_ctc"]
     del mf_ctc["block_id"]  # we don't write this
     assert isinstance(py_ctc["decoupler"], sparse.csc_matrix)
@@ -916,13 +914,8 @@ def test_cross_talk(tmp_path):
     with pytest.warns(RuntimeWarning, match="Not all cross-talk channels"):
         maxwell_filter(raw_missing, cross_talk=ctc_fname)
     # MEG channels not in cross-talk
-    pytest.raises(
-        RuntimeError,
-        maxwell_filter,
-        raw_ctf,
-        origin=(0.0, 0.0, 0.04),
-        cross_talk=ctc_fname,
-    )
+    with pytest.raises(RuntimeError, match="Missing MEG channels"):
+        maxwell_filter(raw_ctf, origin=(0.0, 0.0, 0.04), cross_talk=ctc_fname)
 
 
 @testing.requires_testing_data
@@ -970,10 +963,10 @@ def test_head_translation():
         read_info(sample_fname)["dev_head_t"]["trans"],
     )
     # Degenerate cases
-    pytest.raises(
-        RuntimeError, maxwell_filter, raw, destination=mf_head_origin, coord_frame="meg"
-    )
-    pytest.raises(ValueError, maxwell_filter, raw, destination=[0.0] * 4)
+    with pytest.raises(RuntimeError, match=".* can only be set .* head .*"):
+        maxwell_filter(raw, destination=mf_head_origin, coord_frame="meg")
+    with pytest.raises(ValueError, match="destination must be"):
+        maxwell_filter(raw, destination=[0.0] * 4)
 
 
 # TODO: Eventually add simulation tests mirroring Taulu's original paper
@@ -1395,7 +1388,7 @@ def test_all():
     coord_frames = ("head", "head", "meg", "head")
     ctcs = (ctc_fname, ctc_fname, ctc_fname, ctc_mgh_fname)
     mins = (3.5, 3.5, 1.2, 0.9)
-    meds = (10.8, 10.4, 3.2, 6.0)
+    meds = (10.8, 10.2, 3.2, 5.9)
     st_durs = (1.0, 1.0, 1.0, None)
     destinations = (None, sample_fname, None, None)
     origins = (mf_head_origin, mf_head_origin, mf_meg_origin, mf_head_origin)
@@ -1436,7 +1429,7 @@ def test_triux():
     sss_py = maxwell_filter(
         raw, coord_frame="meg", regularize=None, calibration=tri_cal_fname
     )
-    assert_meg_snr(sss_py, read_crop(tri_sss_cal_fname), 22, 200)
+    assert_meg_snr(sss_py, read_crop(tri_sss_cal_fname), 5, 100)
     # ctc+cal
     sss_py = maxwell_filter(
         raw,
@@ -1445,7 +1438,7 @@ def test_triux():
         calibration=tri_cal_fname,
         cross_talk=tri_ctc_fname,
     )
-    assert_meg_snr(sss_py, read_crop(tri_sss_ctc_cal_fname), 28, 200)
+    assert_meg_snr(sss_py, read_crop(tri_sss_ctc_cal_fname), 5, 100)
     # regularization
     sss_py = maxwell_filter(raw, coord_frame="meg", regularize="in")
     sss_mf = read_crop(tri_sss_reg_fname)