diff --git a/doc/changes/dev/13619.bugfix.rst b/doc/changes/dev/13619.bugfix.rst new file mode 100644 index 00000000000..f5f90275367 --- /dev/null +++ b/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`_. diff --git a/mne/forward/_make_forward.py b/mne/forward/_make_forward.py index 9c2f2552cc9..98859fbcfee 100644 --- a/mne/forward/_make_forward.py +++ b/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, diff --git a/mne/forward/forward.py b/mne/forward/forward.py index 83aa4b3a9b8..e27422d9e05 100644 --- a/mne/forward/forward.py +++ b/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: diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index 9e94f0cbd44..6faa03dac49 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/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): diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 438003c16ee..9dae609f376 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -25,6 +25,7 @@ combine_evoked, compute_raw_covariance, convert_forward_solution, + create_info, make_ad_hoc_cov, make_forward_solution, make_sphere_model, @@ -34,7 +35,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 @@ -1712,3 +1715,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]) diff --git a/mne/simulation/raw.py b/mne/simulation/raw.py index 518e0ffe451..9babcc57afc 100644 --- a/mne/simulation/raw.py +++ b/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, diff --git a/mne/simulation/tests/test_raw.py b/mne/simulation/tests/test_raw.py index 4b1d514e8d7..9babdb7a777 100644 --- a/mne/simulation/tests/test_raw.py +++ b/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():