diff --git a/doc/changes/dev/13068.bugfix.rst b/doc/changes/dev/13068.bugfix.rst new file mode 100644 index 00000000000..91795f1e730 --- /dev/null +++ b/doc/changes/dev/13068.bugfix.rst @@ -0,0 +1 @@ +Fixed ICA getting sources for concatenated raw instances, by :newcontrib:`Beige Jin`. \ No newline at end of file diff --git a/doc/changes/dev/13307.newfeature.rst b/doc/changes/dev/13307.newfeature.rst new file mode 100644 index 00000000000..5917c49531f --- /dev/null +++ b/doc/changes/dev/13307.newfeature.rst @@ -0,0 +1 @@ +Added ``on_inside="raise"`` parameter to :func:`mne.make_forward_solution` and :func:`mne.make_forward_dipole` to control behavior when MEG sensors are inside the outer skin surface. This is useful for forward solutions that are computed with sensors just inside the outer skin surface (e.g., with some OPM coregistrations), by `Eric Larson`_. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index d014a50cab0..2d362aeac66 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -33,6 +33,7 @@ .. _Ashley Drew: https://github.com/ashdrew .. _Asish Panda: https://github.com/kaichogami .. _Austin Hurst: https://github.com/a-hurst +.. _Beige Jin: https://github.com/BeiGeJin .. _Ben Beasley: https://github.com/musicinmybrain .. _Bradley Voytek: https://github.com/voytek .. _Britta Westner: https://britta-wstnr.github.io diff --git a/mne/forward/_make_forward.py b/mne/forward/_make_forward.py index d9447078d29..e1264bd2fe6 100644 --- a/mne/forward/_make_forward.py +++ b/mne/forward/_make_forward.py @@ -37,7 +37,15 @@ apply_trans, invert_transform, ) -from ..utils import _check_fname, _pl, _validate_type, logger, verbose, warn +from ..utils import ( + _check_fname, + _on_missing, + _pl, + _validate_type, + logger, + verbose, + warn, +) from ._compute_forward import ( _compute_forwards, _compute_forwards_meeg, @@ -437,6 +445,7 @@ def _prepare_for_forward( bem, mindist, n_jobs, + *, bem_extra="", trans="", info_extra="", @@ -444,6 +453,7 @@ def _prepare_for_forward( eeg=True, ignore_ref=False, allow_bem_none=False, + on_inside="raise", verbose=None, ): """Prepare for forward computation. @@ -567,11 +577,12 @@ def check_inside_head(x): meg_loc = apply_trans(invert_transform(mri_head_t), meg_loc) n_inside = check_inside_head(meg_loc).sum() if n_inside: - raise RuntimeError( + msg = ( f"Found {n_inside} MEG sensor{_pl(n_inside)} inside the " f"{check_surface}, perhaps coordinate frames and/or " - "coregistration must be incorrect" + "coregistration are incorrect" ) + _on_missing(on_inside, msg, name="on_inside", error_klass=RuntimeError) if len(src): rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) @@ -610,6 +621,7 @@ def make_forward_solution( mindist=0.0, ignore_ref=False, n_jobs=None, + on_inside="raise", verbose=None, ): """Calculate a forward solution for a subject. @@ -640,6 +652,13 @@ def make_forward_solution( option should be True for KIT files, since forward computation with reference channels is not currently supported. %(n_jobs)s + on_inside : 'raise' | 'warn' | 'ignore' + What to do if MEG sensors are inside the outer skin surface. If 'raise' + (default), an error is raised. If 'warn' or 'ignore', the forward + solution is computed anyway and a warning is or isn't emitted, + respectively. + + .. versionadded:: 1.10 %(verbose)s Returns @@ -710,12 +729,13 @@ def make_forward_solution( bem, mindist, n_jobs, - bem_extra, - trans, - info_extra, - meg, - eeg, - ignore_ref, + bem_extra=bem_extra, + trans=trans, + info_extra=info_extra, + meg=meg, + eeg=eeg, + ignore_ref=ignore_ref, + on_inside=on_inside, ) del (src, mri_head_t, trans, info_extra, bem_extra, mindist, meg, eeg, ignore_ref) @@ -741,7 +761,9 @@ def make_forward_solution( @verbose -def make_forward_dipole(dipole, bem, info, trans=None, n_jobs=None, *, verbose=None): +def make_forward_dipole( + dipole, bem, info, trans=None, n_jobs=None, *, on_inside="raise", verbose=None +): """Convert dipole object to source estimate and calculate forward operator. The instance of Dipole is converted to a discrete source space, @@ -767,6 +789,13 @@ def make_forward_dipole(dipole, bem, info, trans=None, n_jobs=None, *, verbose=N The head<->MRI transform filename. Must be provided unless BEM is a sphere model. %(n_jobs)s + on_inside : 'raise' | 'warn' | 'ignore' + What to do if MEG sensors are inside the outer skin surface. If 'raise' + (default), an error is raised. If 'warn' or 'ignore', the forward + solution is computed anyway and a warning is or isn't emitted, + respectively. + + .. versionadded:: 1.10 %(verbose)s Returns @@ -805,7 +834,9 @@ def make_forward_dipole(dipole, bem, info, trans=None, n_jobs=None, *, verbose=N # Forward operator created for channels in info (use pick_info to restrict) # Use defaults for most params, including min_dist - fwd = make_forward_solution(info, trans, src, bem, n_jobs=n_jobs, verbose=verbose) + fwd = make_forward_solution( + info, trans, src, bem, n_jobs=n_jobs, on_inside=on_inside, verbose=verbose + ) # Convert from free orientations to fixed (in-place) convert_forward_solution( fwd, surf_ori=False, force_fixed=True, copy=False, use_cps=False, verbose=None diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index 268dd263af9..b584202465c 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -894,8 +894,11 @@ def test_sensors_inside_bem(): trans["trans"][2, 3] = 0.03 sphere_noshell = make_sphere_model((0.0, 0.0, 0.0), None) sphere = make_sphere_model((0.0, 0.0, 0.0), 1.01) - with pytest.raises(RuntimeError, match=".* 15 MEG.*inside the scalp.*"): - make_forward_solution(info, trans, fname_src, fname_bem) + with pytest.warns(RuntimeWarning, match=".* 15 MEG.*inside the scalp.*"): + fwd = make_forward_solution(info, trans, fname_src, fname_bem, on_inside="warn") + assert fwd["nsource"] == 516 + assert fwd["nchan"] == 42 + assert np.isfinite(fwd["sol"]["data"]).all() make_forward_solution(info, trans, fname_src, fname_bem_meg) # okay make_forward_solution(info, trans, fname_src, sphere_noshell) # okay with pytest.raises(RuntimeError, match=".* 42 MEG.*outermost sphere sh.*"): diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index 6805007927e..0876d40911a 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -1295,8 +1295,10 @@ def _sources_as_raw(self, raw, add_channels, start, stop): picks = pick_channels(raw.ch_names, add_channels) data_ = np.concatenate([data_, raw.get_data(picks, start=start, stop=stop)]) out._data = data_ - out._first_samps = [out.first_samp] - out._last_samps = [out.last_samp] + out_first_samp = out.first_samp + out_last_samp = out.last_samp + out._first_samps = [out_first_samp] + out._last_samps = [out_last_samp] out.filenames = [None] out.preload = True out._projector = None diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py index ab8d401e019..b56e214b850 100644 --- a/mne/preprocessing/tests/test_ica.py +++ b/mne/preprocessing/tests/test_ica.py @@ -25,6 +25,7 @@ EpochsArray, EvokedArray, Info, + concatenate_raws, create_info, make_ad_hoc_cov, pick_channels_regexp, @@ -1727,3 +1728,21 @@ def test_ica_ch_types(ch_type): for inst in [raw, epochs, evoked]: ica.apply(inst) ica.get_sources(inst) + + +@testing.requires_testing_data +def test_ica_get_sources_concatenated(): + """Test ICA get_sources method with concatenated raws.""" + # load data + raw = read_raw_fif(raw_fname).crop(0, 3).load_data() # raw has 3 seconds of data + # create concatenated raw instances + raw_concat = concatenate_raws( + [raw.copy(), raw.copy()] + ) # raw_concat has 6 seconds of data + # do ICA + ica = ICA(n_components=2, max_iter=2) + with _record_warnings(), pytest.warns(UserWarning, match="did not converge"): + ica.fit(raw_concat) + # get sources + raw_sources = ica.get_sources(raw_concat) # but this only has 3 seconds of data + assert raw_concat.n_times == raw_sources.n_times # this will fail diff --git a/mne/simulation/raw.py b/mne/simulation/raw.py index 94f186992e1..518e0ffe451 100644 --- a/mne/simulation/raw.py +++ b/mne/simulation/raw.py @@ -46,6 +46,7 @@ _check_preload, _pl, _validate_type, + _verbose_safe_false, check_random_state, logger, verbose, @@ -793,7 +794,14 @@ def _iter_forward_solutions( info.update(projs=[], bads=[]) # Ensure no 'projs' or 'bads' mri_head_t, trans = _get_trans(trans) sensors, rr, info, update_kwargs, bem = _prepare_for_forward( - src, mri_head_t, info, bem, mindist, n_jobs, allow_bem_none=True, verbose=False + src, + mri_head_t, + info, + bem, + mindist, + n_jobs, + allow_bem_none=True, + verbose=_verbose_safe_false(), ) del (src, mindist)