BUG: Fix bug with fwd.info["dev_head_t"] setting (#13619)

Author: larsoner

doc/changes/dev/13619.bugfix.rst

@@ -0,0 +1 @@
+Fix bug with :func:`mne.make_forward_solution` when no MEG channels are present where ``dev_head_t`` was errantly set, by `Eric Larson`_.

mne/forward/_make_forward.py

@@ -499,8 +499,7 @@ def _prepare_for_forward(
     kwargs_fwd_info = dict(
         chs=info["chs"],
         comps=info["comps"],
-        # The forward-writing code always wants a dev_head_t, so give an identity one
-        dev_head_t=info["dev_head_t"] or Transform("meg", "head"),
+        dev_head_t=info["dev_head_t"],
         mri_file=info_trans,
         mri_id=mri_id,
         meas_file=info_extra,

mne/forward/forward.py

@@ -441,14 +441,13 @@ def _read_forward_meas_info(tree, fid):
     # Get the MEG device <-> head coordinate transformation
     tag = find_tag(fid, parent_meg, FIFF.FIFF_COORD_TRANS)
     if tag is None:
-        raise ValueError("MEG/head coordinate transformation not found")
-    cand = tag.data
-    if cand["from"] == coord_device and cand["to"] == coord_head:
-        info["dev_head_t"] = cand
-    elif cand["from"] == coord_ctf_head and cand["to"] == coord_head:
-        info["ctf_head_t"] = cand
+        info["dev_head_t"] = None  # EEG-only
     else:
-        raise ValueError("MEG/head coordinate transformation not found")
+        cand = tag.data
+        if cand["from"] == coord_device and cand["to"] == coord_head:
+            info["dev_head_t"] = cand
+        elif cand["from"] == coord_ctf_head and cand["to"] == coord_head:
+            info["ctf_head_t"] = cand
 
     bads = _read_bad_channels(fid, parent_meg, ch_names_mapping=ch_names_mapping)
     # clean up our bad list, old versions could have non-existent bads
@@ -1107,9 +1106,8 @@ def write_forward_meas_info(fid, info):
         write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, info["meas_id"])
     # get transformation from CTF and DEVICE to HEAD coordinate frame
     meg_head_t = info.get("dev_head_t", info.get("ctf_head_t"))
-    if meg_head_t is None:
-        raise ValueError("Head<-->sensor transform not found")
-    write_coord_trans(fid, meg_head_t)
+    if meg_head_t is not None:
+        write_coord_trans(fid, meg_head_t)
 
     ch_names_mapping = dict()
     if "chs" in info:

mne/forward/tests/test_make_forward.py

@@ -836,11 +836,14 @@ def test_make_forward_no_meg(tmp_path):
     trans = None
     montage = make_standard_montage("standard_1020")
     info = create_info(["Cz"], 1000.0, "eeg").set_montage(montage)
+    assert info["dev_head_t"] is None  # gh-13604
     fwd = make_forward_solution(info, trans, src, bem)
+    assert fwd["info"]["dev_head_t"] is None
     fname = tmp_path / "test-fwd.fif"
     write_forward_solution(fname, fwd)
     fwd_read = read_forward_solution(fname)
     assert_allclose(fwd["sol"]["data"], fwd_read["sol"]["data"])
+    assert fwd_read["info"]["dev_head_t"] is None
 
 
 def test_use_coil_def(tmp_path):

mne/minimum_norm/tests/test_inverse.py

@@ -26,6 +26,7 @@
     compute_rank,
     compute_raw_covariance,
     convert_forward_solution,
+    create_info,
     make_ad_hoc_cov,
     make_forward_solution,
     make_sphere_model,
@@ -35,7 +36,9 @@
     read_cov,
     read_evokeds,
     read_forward_solution,
+    setup_volume_source_space,
 )
+from mne.channels import make_standard_montage
 from mne.datasets import testing
 from mne.epochs import Epochs, EpochsArray, make_fixed_length_epochs
 from mne.event import read_events
@@ -1726,3 +1729,38 @@ def test_allow_mixed_source_spaces(mixed_fwd_cov_evoked):
     inv_op = make_inverse_operator(evoked.info, fwd, cov, loose=0.0)
     stc = apply_inverse(evoked, inv_op, lambda2=1.0 / 9.0)  # normal
     assert (stc.data < 0).any()
+
+
+@testing.requires_testing_data
+def test_make_inverse_no_meg(tmp_path):
+    """Test that we can make and I/O forward solution with no MEG channels."""
+    pos = dict(rr=[[0.05, 0, 0]], nn=[[0, 0, 1.0]])
+    src = setup_volume_source_space(pos=pos)
+    bem = make_sphere_model()
+    trans = None
+    montage = make_standard_montage("standard_1020")
+    info = create_info(["Cz", "Fz"], 1000.0, "eeg").set_montage(montage)
+    data = np.random.default_rng(0).standard_normal((len(info["ch_names"]), 10))
+    evoked = EvokedArray(data, info)
+    evoked.set_eeg_reference(projection=True)
+    del info, data
+    fwd = make_forward_solution(evoked.info, trans, src, bem)
+    assert fwd["info"]["dev_head_t"] is None
+    cov = make_ad_hoc_cov(evoked.info)
+    # Ensure we maintain backward-compatible behavior with when we *did* save
+    # the forward with an identity dev_head_t
+    fname = tmp_path / "test-inv.fif"
+    invs = list()
+    for dev_head_t in (None, mne.Transform("meg", "head")):
+        fwd["info"]["dev_head_t"] = dev_head_t
+        inv = make_inverse_operator(evoked.info, fwd, cov)
+        assert inv["info"]["dev_head_t"] == dev_head_t
+        stc = apply_inverse(evoked, inv, lambda2=1.0 / 9.0)
+        assert stc.data.shape == (1, 10)
+        write_inverse_operator(fname, inv, overwrite=True)
+        inv_read = read_inverse_operator(fname)
+        assert inv_read["info"]["dev_head_t"] == dev_head_t
+        invs.append(inv)
+    assert invs[1]["info"]["dev_head_t"] is not None
+    invs[1]["info"]["dev_head_t"] = None  # for comparison
+    _compare(invs[0], invs[1])

mne/simulation/raw.py

@@ -41,7 +41,7 @@
     setup_volume_source_space,
 )
 from ..surface import _CheckInside
-from ..transforms import _get_trans, transform_surface_to
+from ..transforms import Transform, _get_trans, transform_surface_to
 from ..utils import (
     _check_preload,
     _pl,
@@ -131,15 +131,16 @@ def _check_head_pos(head_pos, info, first_samp, times=None):
                 f"s), found {bad.sum()}/{len(bad)} bad values (is this a split "
                 "file?)"
             )
+    dev_head_t = info["dev_head_t"] or Transform("meg", "head")  # deal with None
     # If it starts close to zero, make it zero (else unique(offset) fails)
     if len(ts) > 0 and ts[0] < (0.5 / info["sfreq"]):
         ts[0] = 0.0
     # If it doesn't start at zero, insert one at t=0
     elif len(ts) == 0 or ts[0] > 0:
         ts = np.r_[[0.0], ts]
-        dev_head_ts.insert(0, info["dev_head_t"]["trans"])
+        dev_head_ts.insert(0, dev_head_t["trans"])
     dev_head_ts = [
-        {"trans": d, "to": info["dev_head_t"]["to"], "from": info["dev_head_t"]["from"]}
+        {"trans": d, "to": dev_head_t["to"], "from": dev_head_t["from"]}
         for d in dev_head_ts
     ]
     offsets = np.round(ts * info["sfreq"]).astype(int)
@@ -290,7 +291,8 @@ def simulate_raw(
                 "If forward is not None then trans, src, bem, "
                 "and head_pos must all be None"
             )
-        if not np.allclose(
+        ch_types = info.get_channel_types(unique=True)
+        if ("mag" in ch_types or "grad" in ch_types) and not np.allclose(
             forward["info"]["dev_head_t"]["trans"],
             info["dev_head_t"]["trans"],
             atol=1e-6,

mne/simulation/tests/test_raw.py

@@ -54,7 +54,7 @@
 from mne.source_space._source_space import _compare_source_spaces
 from mne.surface import _get_ico_surface
 from mne.tests.test_chpi import _assert_quats
-from mne.transforms import _affine_to_quat
+from mne.transforms import Transform, _affine_to_quat
 from mne.utils import catch_logging
 
 raw_fname_short = Path(__file__).parents[2] / "io" / "tests" / "data" / "test_raw.fif"
@@ -501,6 +501,25 @@ def test_simulate_round_trip(raw_data):
         simulate_raw(raw.info, stc, None, None, None, forward=fwd)
 
 
+def test_simulate_eeg_only(raw_data):
+    """Test simulate_raw with EEG data only."""
+    # gh-13604
+    raw, src, stc, trans, sphere = raw_data
+    raw = raw.pick("eeg")
+    assert raw.info["dev_head_t"] is not None
+    with raw.info._unlock():
+        raw.info["dev_head_t"] = None
+    fwd = make_forward_solution(raw.info, trans, src, sphere)
+    assert fwd["info"]["dev_head_t"] is None
+    # forward with identity dev_head_t (old style) but EEG-only data
+    raw_sim = simulate_raw(raw.info, stc, None, None, None, forward=fwd)
+    fwd["info"]["dev_head_t"] = Transform("meg", "head")
+    raw_sim_2 = simulate_raw(raw.info, stc, None, None, None, forward=fwd)
+    assert raw_sim.info["dev_head_t"] is None
+    assert raw_sim_2.info["dev_head_t"] is None
+    assert_allclose(raw_sim[:][0], raw_sim_2[:][0])
+
+
 @pytest.mark.slowtest
 @testing.requires_testing_data
 def test_simulate_raw_chpi():