From 909fa85669bbf719d7f506f76b7c6faea8e9c5f1 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 3 Sep 2025 17:10:32 -0400 Subject: [PATCH 1/8] ENH: Speed up forward computations for iterative fitting --- doc/changes/dev/bugfix.rst | 1 + mne/dipole.py | 2 +- mne/forward/_compute_forward.py | 4 +- mne/forward/_make_forward.py | 108 ++++++++++++++++++++++--- mne/forward/tests/test_make_forward.py | 23 ++++++ mne/source_space/_source_space.py | 49 ++++++----- mne/surface.py | 39 ++++++++- 7 files changed, 189 insertions(+), 37 deletions(-) create mode 100644 doc/changes/dev/bugfix.rst diff --git a/doc/changes/dev/bugfix.rst b/doc/changes/dev/bugfix.rst new file mode 100644 index 00000000000..150d0bbfe3f --- /dev/null +++ b/doc/changes/dev/bugfix.rst @@ -0,0 +1 @@ +Fix bug with :func:`mne.make_forward_solution` where sources were not checked to make sure they're inside the inner skull for spherical BEMs, by `Eric Larson`_. diff --git a/mne/dipole.py b/mne/dipole.py index 6978fe77936..4fb1d9aa992 100644 --- a/mne/dipole.py +++ b/mne/dipole.py @@ -1720,7 +1720,7 @@ def fit_dipole( # C code computes guesses w/sphere model for speed, don't bother here fwd_data = _prep_field_computation( - guess_src["rr"], sensors=sensors, bem=bem, n_jobs=n_jobs, verbose=safe_false + sensors=sensors, bem=bem, n_jobs=n_jobs, verbose=safe_false ) fwd_data["inner_skull"] = inner_skull guess_fwd, guess_fwd_orig, guess_fwd_scales = _dipole_forwards( diff --git a/mne/forward/_compute_forward.py b/mne/forward/_compute_forward.py index 12358a0cbe7..8ba15def389 100644 --- a/mne/forward/_compute_forward.py +++ b/mne/forward/_compute_forward.py @@ -709,7 +709,7 @@ def _compute_mdfv(rrs, rmags, cosmags, ws, bins, too_close): @verbose -def _prep_field_computation(rr, *, sensors, bem, n_jobs, verbose=None): +def _prep_field_computation(*, sensors, bem, n_jobs, verbose=None): """Precompute and store some things that are used for both MEG and EEG. Calculation includes multiplication factors, coordinate transforms, @@ -840,7 +840,7 @@ def _compute_forwards(rr, *, bem, sensors, n_jobs, verbose=None): # This modifies "sensors" in place, so let's copy it in case the calling # function needs to reuse it (e.g., in simulate_raw.py) sensors = deepcopy(sensors) - fwd_data = _prep_field_computation(rr, sensors=sensors, bem=bem, n_jobs=n_jobs) + fwd_data = _prep_field_computation(sensors=sensors, bem=bem, n_jobs=n_jobs) Bs = _compute_forwards_meeg( rr, sensors=sensors, fwd_data=fwd_data, n_jobs=n_jobs ) diff --git a/mne/forward/_make_forward.py b/mne/forward/_make_forward.py index 421cb921c1d..3fa08b87183 100644 --- a/mne/forward/_make_forward.py +++ b/mne/forward/_make_forward.py @@ -26,7 +26,7 @@ _filter_source_spaces, _make_discrete_source_space, ) -from ..surface import _CheckInside, _normalize_vectors +from ..surface import _CheckInside, _CheckInsideSphere, _normalize_vectors from ..transforms import ( Transform, _coord_frame_name, @@ -38,7 +38,11 @@ transform_surface_to, ) from ..utils import _check_fname, _pl, _validate_type, logger, verbose, warn -from ._compute_forward import _compute_forwards +from ._compute_forward import ( + _compute_forwards, + _compute_forwards_meeg, + _prep_field_computation, +) from .forward import _FWD_ORDER, Forward, _merge_fwds, convert_forward_solution _accuracy_dict = dict( @@ -533,35 +537,38 @@ def _prepare_for_forward( # Circumvent numerical problems by excluding points too close to the skull, # and check that sensors are not inside any BEM surface if bem is not None: + kwargs = dict(limit=mindist, mri_head_t=mri_head_t, src=src) if not bem["is_sphere"]: check_surface = "inner skull surface" - inner_skull = _bem_find_surface(bem, "inner_skull") - check_inside = _filter_source_spaces( - inner_skull, mindist, mri_head_t, src, n_jobs - ) + check_inside_brain = _CheckInside(_bem_find_surface(bem, "inner_skull")) + _filter_source_spaces(check_inside_brain, n_jobs=n_jobs, **kwargs) logger.info("") if len(bem["surfs"]) == 3: check_surface = "scalp surface" - check_inside = _CheckInside(_bem_find_surface(bem, "head")) + check_inside_head = _CheckInside(_bem_find_surface(bem, "head")) else: check_surface = "outermost sphere shell" + check_inside_brain = _CheckInsideSphere(bem) + _filter_source_spaces(check_inside_brain, **kwargs) if len(bem["layers"]) == 0: - def check_inside(x): + def check_inside_head(x): return np.zeros(len(x), bool) else: - def check_inside(x): + def check_inside_head(x): + # MEG sensors are in MRI coords at this point, so need to + # transform BEM center to MRI coords, too r0 = apply_trans(invert_transform(mri_head_t), bem["r0"]) - return np.linalg.norm(x - r0, axis=1) < bem["layers"][-1]["rad"] + return np.linalg.norm(x - r0, axis=1) < bem.radius if "meg" in sensors: meg_loc = apply_trans( invert_transform(mri_head_t), np.array([coil["r0"] for coil in sensors["meg"]["defs"]]), ) - n_inside = check_inside(meg_loc).sum() + n_inside = check_inside_head(meg_loc).sum() if n_inside: raise RuntimeError( f"Found {n_inside} MEG sensor{_pl(n_inside)} inside the " @@ -934,3 +941,82 @@ def use_coil_def(fname): yield finally: _extra_coil_def_fname = None + + +# Class optimized for incremental fitting using the same sensors and BEM at different +# locations (e.g., for Xfit-like iterative forward solutions) +class _ForwardModeler: + @verbose + def __init__( + self, + info, + trans, + bem, + *, + mindist=0.0, + n_jobs=1, + verbose=None, + ): + self.mri_head_t, _ = _get_trans(trans) + self.mindist = mindist + self.n_jobs = n_jobs + # TODO: Make `src` optional in _prepare_for_forward + import mne + + src = mne.setup_volume_source_space( + "", pos=dict(rr=np.zeros((1, 3)), nn=np.array([[0, 0, 1]])), verbose="error" + ) + self.sensors, _, _, _, self.bem = _prepare_for_forward( + src, + self.mri_head_t, + info, + bem, + mindist, + n_jobs, + bem_extra="", + trans="", + info_extra="", + meg=True, + eeg=True, + ignore_ref=False, + ) + self.fwd_data = _prep_field_computation( + sensors=self.sensors, + bem=self.bem, + n_jobs=self.n_jobs, + ) + if self.bem["is_sphere"]: + self.check_inside = _CheckInsideSphere(self.bem) + else: + self.check_inside = _CheckInside(_bem_find_surface(self.bem, "inner_skull")) + + def compute(self, src): + src = _ensure_src(src).copy() + for s in src: + transform_surface_to(s, "head", self.mri_head_t) + + kwargs = dict(limit=self.mindist, mri_head_t=self.mri_head_t, src=src) + _filter_source_spaces(self.check_inside, n_jobs=self.n_jobs, **kwargs) + rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) + if len(rr) < 1: + raise RuntimeError( + "No points left in source space after excluding " + "points close to inner skull." + ) + + sensors = deepcopy(self.sensors) + fwd_data = deepcopy(self.fwd_data) + fwds = _compute_forwards_meeg( + rr, + sensors=sensors, + fwd_data=fwd_data, + n_jobs=self.n_jobs, + ) + fwds = { + key: _to_forward_dict(fwds[key], sensors[key]["ch_names"]) + for key in _FWD_ORDER + if key in fwds + } + fwd = _merge_fwds(fwds, verbose=False) + del fwds + return fwd diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index d408df3b551..9e70b08a013 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -602,6 +602,29 @@ def test_make_forward_solution_sphere(tmp_path, fname_src_small): assert fwd["mri_head_t"]["trans"][0, 3] == -0.05 +@testing.requires_testing_data +def test_make_forward_sphere_exclude(): + """Test that points are excluded that are outside BEM sphere inner layer.""" + bem = make_sphere_model(r0=(0.0, 0.0, 0.04), head_radius=0.1) + src = setup_volume_source_space( + "sample", + pos=10.0, + sphere=(0.0, 0.0, 0.04, 0.1), + subjects_dir=subjects_dir, + exclude=10, + verbose=True, + ) + assert src[0]["nuse"] == 3694 + trans = Transform("mri", "head") # identity for simplicity + raw = read_raw_fif(fname_raw) + raw.pick(raw.ch_names[:1]) + fwd = make_forward_solution(raw.info, trans, src, bem, verbose=True) + assert fwd["nsource"] == 3694 + bem_small = make_sphere_model(r0=(0.0, 0.0, 0.04), head_radius=0.09) + fwd_small = make_forward_solution(raw.info, trans, src, bem_small, verbose=True) + assert fwd_small["nsource"] < 3694 + + @pytest.mark.slowtest @testing.requires_testing_data def test_forward_mixed_source_space(tmp_path): diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 7919b095f28..7e2e27d1961 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -45,6 +45,7 @@ from ..parallel import parallel_func from ..surface import ( _CheckInside, + _CheckInsideSphere, _compute_nearest, _create_surf_spacing, _get_ico_surface, @@ -1907,7 +1908,9 @@ def setup_volume_source_space( surf["rr"] *= 1e-3 # must be converted to meters else: # Load an icosahedron and use that as the surface logger.info("Setting up the sphere...") - surf = dict(R=sphere[3], r0=sphere[:3]) + surf = ConductorModel( + layers=[dict(rad=sphere[3])], r0=sphere[:3], is_sphere=True + ) # Make the grid of sources in MRI space sp = _make_volume_source_space( surf, @@ -2063,10 +2066,10 @@ def _make_volume_source_space( cm = np.mean(surf["rr"], axis=0) # center of mass maxdist = np.linalg.norm(surf["rr"] - cm, axis=1).max() else: - mins = surf["r0"] - surf["R"] - maxs = surf["r0"] + surf["R"] + mins = surf["r0"] - surf.radius + maxs = surf["r0"] + surf.radius cm = surf["r0"].copy() - maxdist = surf["R"] + maxdist = surf.radius # Define the sphere which fits the surface logger.info( @@ -2136,18 +2139,13 @@ def _make_volume_source_space( 1000 * exclude, 1000 * maxdist, ) + kwargs = dict(limit=mindist, mri_head_t=None, src=[sp]) + assert sp["coord_frame"] == FIFF.FIFFV_COORD_MRI if "rr" in surf: - _filter_source_spaces(surf, mindist, None, [sp], n_jobs) - else: # sphere - vertno = np.where(sp["inuse"])[0] - bads = ( - np.linalg.norm(sp["rr"][vertno] - surf["r0"], axis=-1) - >= surf["R"] - mindist / 1000.0 - ) - sp["nuse"] -= bads.sum() - sp["inuse"][vertno[bads]] = False - sp["vertno"] = np.where(sp["inuse"])[0] - del vertno + assert surf["coord_frame"] == FIFF.FIFFV_COORD_MRI + else: + assert surf["is_sphere"] + _filter_source_spaces(surf, n_jobs=n_jobs, **kwargs) del surf logger.info( "%d sources remaining after excluding the sources outside " @@ -2534,7 +2532,9 @@ def _grid_interp_jit(from_shape, to_shape, trans, order, inuse): @verbose -def _filter_source_spaces(surf, limit, mri_head_t, src, n_jobs=None, verbose=None): +def _filter_source_spaces( + surf_or_check_inside, *, limit, mri_head_t, src, n_jobs=None, verbose=None +): """Remove all source space points closer than a given limit (in mm).""" if src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD and mri_head_t is None: raise RuntimeError( @@ -2558,7 +2558,14 @@ def _filter_source_spaces(surf, limit, mri_head_t, src, n_jobs=None, verbose=Non logger.info(out_str + " (will take a few...)") # fit a sphere to a surf quickly - check_inside = _CheckInside(surf) + if isinstance(surf_or_check_inside, _CheckInside | _CheckInsideSphere): + check_inside = surf_or_check_inside + else: + if "rr" in surf_or_check_inside: + check_inside = _CheckInside(surf_or_check_inside) + else: + check_inside = _CheckInsideSphere(surf_or_check_inside) + del surf_or_check_inside # Check that the source is inside surface (often the inner skull) for s in src: @@ -2568,7 +2575,7 @@ def _filter_source_spaces(surf, limit, mri_head_t, src, n_jobs=None, verbose=Non if s["coord_frame"] == FIFF.FIFFV_COORD_HEAD: r1s = apply_trans(inv_trans["trans"], r1s) - inside = check_inside(r1s, n_jobs) + inside = check_inside(r1s, n_jobs=n_jobs) omit_outside = (~inside).sum() # vectorized nearest using BallTree (or cdist) @@ -2578,16 +2585,14 @@ def _filter_source_spaces(surf, limit, mri_head_t, src, n_jobs=None, verbose=Non idx = np.where(inside)[0] check_r1s = r1s[idx] if check_inside.inner_r is not None: - # ... and those that are at least inner_sphere + limit away + # ... and those that are at least inner_sphere - limit away mask = ( np.linalg.norm(check_r1s - check_inside.cm, axis=-1) >= check_inside.inner_r - limit / 1000.0 ) idx = idx[mask] check_r1s = check_r1s[mask] - dists = _compute_nearest( - surf["rr"], check_r1s, return_dists=True, method="KDTree" - )[1] + dists = check_inside.query(check_r1s)[0] close = dists < limit / 1000.0 omit_limit = np.sum(close) inside[idx[close]] = False diff --git a/mne/surface.py b/mne/surface.py index 1a2b05efe81..23367974e86 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -779,7 +779,7 @@ def _init_pyvista(self): self.pdata = _surface_to_polydata(self.surf).clean() @verbose - def __call__(self, rr, n_jobs=None, verbose=None): + def __call__(self, rr, *, n_jobs=None, verbose=None): n_orig = len(rr) logger.info( f"Checking surface interior status for {n_orig} point{_pl(n_orig, ' ')}..." @@ -794,6 +794,14 @@ def __call__(self, rr, n_jobs=None, verbose=None): logger.info(f"Interior check completed in {(time.time() - t0) * 1000:0.1f} ms") return inside + def query(self, rr): + """Get the distance to the nearest point.""" + if not hasattr(self, "_tree"): # compute on the fly only when needed + from scipy.spatial import KDTree + + self._tree = KDTree(self.surf["rr"]) + return self._tree.query(rr) + def _call_pyvista(self, rr): pdata = _surface_to_polydata(dict(rr=rr)) out = pdata.select_enclosed_points(self.pdata, check_surface=False) @@ -855,6 +863,35 @@ def _call_old(self, rr, n_jobs): return inside +class _CheckInsideSphere: + def __init__(self, sphere): + from .bem import ConductorModel + + assert isinstance(sphere, ConductorModel) and sphere["is_sphere"] + self.cm = sphere["r0"] + self.inner_r = sphere.radius # float or None + + # No need for verbose dec here because no MNE code is called that would log + def __call__(self, rr, *, n_jobs=None, verbose=None): + assert isinstance(rr, np.ndarray), type(rr) + assert rr.ndim == 2 and rr.shape[1] == 3 + if self.inner_r is None: + return np.ones(rr.shape[0], bool) + else: + return np.linalg.norm(rr - self.cm, axis=-1) <= self.inner_r + + def query(self, rr): + """Return the distance to the sphere surface for each point.""" + assert isinstance(rr, np.ndarray), type(rr) + assert rr.ndim == 2 and rr.shape[1] == 3, rr.shape + idx = np.zeros(rr.shape[0], int) + if self.inner_r is None: + dists = np.full(rr.shape[0], np.inf) + else: + dists = np.abs(np.linalg.norm(rr - self.cm, axis=-1) - self.inner_r) + return dists, idx + + ############################################################################### # Handle freesurfer From db811116eacefacda74bd2be893c5f6936f274b9 Mon Sep 17 00:00:00 2001 From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com> Date: Wed, 3 Sep 2025 21:21:59 +0000 Subject: [PATCH 2/8] [autofix.ci] apply automated fixes --- doc/changes/dev/{bugfix.rst => 13407.bugfix.rst} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename doc/changes/dev/{bugfix.rst => 13407.bugfix.rst} (100%) diff --git a/doc/changes/dev/bugfix.rst b/doc/changes/dev/13407.bugfix.rst similarity index 100% rename from doc/changes/dev/bugfix.rst rename to doc/changes/dev/13407.bugfix.rst From a1ded17f914ae787f086e5b8032a79586d5643f4 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 4 Sep 2025 11:28:38 -0400 Subject: [PATCH 3/8] FIX: Fixes --- mne/dipole.py | 20 ++++--- mne/forward/_make_forward.py | 62 ++++++++++---------- mne/forward/tests/test_make_forward.py | 79 ++++++++++++++++++++------ mne/source_space/_source_space.py | 4 +- mne/surface.py | 29 +++++++--- mne/tests/test_dipole.py | 13 ++++- 6 files changed, 137 insertions(+), 70 deletions(-) diff --git a/mne/dipole.py b/mne/dipole.py index 4fb1d9aa992..67ff5cf65c4 100644 --- a/mne/dipole.py +++ b/mne/dipole.py @@ -17,7 +17,7 @@ from ._fiff.pick import pick_types from ._fiff.proj import _needs_eeg_average_ref_proj, make_projector from ._freesurfer import _get_aseg, head_to_mni, head_to_mri, read_freesurfer_lut -from .bem import _bem_find_surface, _bem_surf_name, _fit_sphere +from .bem import ConductorModel, _bem_find_surface, _bem_surf_name, _fit_sphere from .cov import _ensure_cov, compute_whitener from .evoked import _aspect_rev, _read_evoked, _write_evokeds from .fixes import _safe_svd @@ -960,9 +960,8 @@ def _make_guesses(surf, grid, exclude, mindist, n_jobs=None, verbose=None): ) else: logger.info( - "Making a spherical guess space with radius {:7.1f} mm...".format( - 1000 * surf["R"] - ) + f"Making a spherical guess space with radius {1000 * surf.radius:7.1f} " + "mm..." ) logger.info("Filtering (grid = %6.f mm)..." % (1000 * grid)) src = _make_volume_source_space( @@ -1317,7 +1316,7 @@ def _fit_dipole( constraint = partial( _sphere_constraint, r0=fwd_data["inner_skull"]["r0"], - R_adj=fwd_data["inner_skull"]["R"] - min_dist_to_inner_skull, + R_adj=fwd_data["inner_skull"].radius - min_dist_to_inner_skull, ) # Find a good starting point (find_best_guess in C) @@ -1600,13 +1599,14 @@ def fit_dipole( raise RuntimeError( "No MEG channels found, but MEG-only sphere model used" ) - R = np.min(np.sqrt(np.sum(R * R, axis=1))) # use dist to sensors - kind = "max_rad" + R = np.min(np.linalg.norm(R, axis=1)) + kind = "min_rad" logger.info( f"Sphere model : origin at ({1000 * r0[0]: 7.2f} {1000 * r0[1]: 7.2f} " f"{1000 * r0[2]: 7.2f}) mm, {kind} = {R:6.1f} mm" ) - inner_skull = dict(R=R, r0=r0) # NB sphere model defined in head frame + # NB sphere model defined in head frame + inner_skull = ConductorModel(layers=[dict(rad=R)], r0=r0, is_sphere=True) del R, r0 # Deal with DipoleFixed cases here @@ -1710,7 +1710,9 @@ def fit_dipole( check = _surface_constraint(pos, inner_skull, min_dist_to_inner_skull) else: check = _sphere_constraint( - pos, inner_skull["r0"], R_adj=inner_skull["R"] - min_dist_to_inner_skull + pos, + inner_skull["r0"], + R_adj=inner_skull.radius - min_dist_to_inner_skull, ) if check <= 0: raise ValueError( diff --git a/mne/forward/_make_forward.py b/mne/forward/_make_forward.py index 3fa08b87183..7342574b704 100644 --- a/mne/forward/_make_forward.py +++ b/mne/forward/_make_forward.py @@ -21,6 +21,7 @@ from ..bem import ConductorModel, _bem_find_surface, read_bem_solution from ..source_estimate import VolSourceEstimate from ..source_space._source_space import ( + SourceSpaces, _complete_vol_src, _ensure_src, _filter_source_spaces, @@ -463,7 +464,7 @@ def _prepare_for_forward( # let's make a copy in case we modify something src = _ensure_src(src).copy() nsource = sum(s["nuse"] for s in src) - if nsource == 0: + if len(src) and nsource == 0: raise RuntimeError( "No sources are active in these source spaces. " '"do_all" option should be used.' @@ -523,9 +524,12 @@ def _prepare_for_forward( # (will either be HEAD or MRI) for s in src: transform_surface_to(s, "head", mri_head_t) - logger.info( - f"Source spaces are now in {_coord_frame_name(s['coord_frame'])} coordinates." - ) + del s + if len(src): + logger.info( + f"Source spaces are now in {_coord_frame_name(src[0]['coord_frame'])} " + "coordinates." + ) # Prepare the BEM model eegnames = sensors.get("eeg", dict()).get("ch_names", []) @@ -541,33 +545,29 @@ def _prepare_for_forward( if not bem["is_sphere"]: check_surface = "inner skull surface" check_inside_brain = _CheckInside(_bem_find_surface(bem, "inner_skull")) - _filter_source_spaces(check_inside_brain, n_jobs=n_jobs, **kwargs) logger.info("") if len(bem["surfs"]) == 3: check_surface = "scalp surface" check_inside_head = _CheckInside(_bem_find_surface(bem, "head")) + else: + check_inside_head = check_inside_brain else: check_surface = "outermost sphere shell" check_inside_brain = _CheckInsideSphere(bem) - _filter_source_spaces(check_inside_brain, **kwargs) - if len(bem["layers"]) == 0: + if bem.radius is not None: + check_inside_head = _CheckInsideSphere(bem, check="outer") + else: def check_inside_head(x): return np.zeros(len(x), bool) - else: - - def check_inside_head(x): - # MEG sensors are in MRI coords at this point, so need to - # transform BEM center to MRI coords, too - r0 = apply_trans(invert_transform(mri_head_t), bem["r0"]) - return np.linalg.norm(x - r0, axis=1) < bem.radius + if len(src): + _filter_source_spaces(check_inside_brain, **kwargs) if "meg" in sensors: - meg_loc = apply_trans( - invert_transform(mri_head_t), - np.array([coil["r0"] for coil in sensors["meg"]["defs"]]), - ) + meg_loc = np.array([coil["r0"] for coil in sensors["meg"]["defs"]]) + if not bem["is_sphere"]: + meg_loc = apply_trans(invert_transform(mri_head_t), meg_loc) n_inside = check_inside_head(meg_loc).sum() if n_inside: raise RuntimeError( @@ -576,12 +576,15 @@ def check_inside_head(x): "coregistration must be incorrect" ) - rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) - if len(rr) < 1: - raise RuntimeError( - "No points left in source space after excluding " - "points close to inner skull." - ) + if len(src): + rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) + if len(rr) < 1: + raise RuntimeError( + "No points left in source space after excluding " + "points close to inner skull." + ) + else: + rr = np.zeros((0, 3)) # deal with free orientations: source_nn = np.tile(np.eye(3), (len(rr), 1)) @@ -943,9 +946,9 @@ def use_coil_def(fname): _extra_coil_def_fname = None -# Class optimized for incremental fitting using the same sensors and BEM at different -# locations (e.g., for Xfit-like iterative forward solutions) class _ForwardModeler: + """Optimized incremental fitting using the same sensors and BEM.""" + @verbose def __init__( self, @@ -960,12 +963,7 @@ def __init__( self.mri_head_t, _ = _get_trans(trans) self.mindist = mindist self.n_jobs = n_jobs - # TODO: Make `src` optional in _prepare_for_forward - import mne - - src = mne.setup_volume_source_space( - "", pos=dict(rr=np.zeros((1, 3)), nn=np.array([[0, 0, 1]])), verbose="error" - ) + src = SourceSpaces([]) self.sensors, _, _, _, self.bem = _prepare_for_forward( src, self.mri_head_t, diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index 9e70b08a013..a82b4cb150b 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -33,7 +33,11 @@ from mne.dipole import Dipole, fit_dipole from mne.forward import Forward, _do_forward_solution, use_coil_def from mne.forward._compute_forward import _magnetic_dipole_field_vec -from mne.forward._make_forward import _create_meg_coils, make_forward_dipole +from mne.forward._make_forward import ( + _create_meg_coils, + _ForwardModeler, + make_forward_dipole, +) from mne.forward.tests.test_forward import assert_forward_allclose from mne.io import read_info, read_raw_bti, read_raw_fif, read_raw_kit from mne.simulation import simulate_evoked @@ -44,7 +48,7 @@ write_source_spaces, ) from mne.surface import _get_ico_surface -from mne.transforms import Transform +from mne.transforms import Transform, apply_trans, invert_transform from mne.utils import ( _record_warnings, catch_logging, @@ -602,27 +606,40 @@ def test_make_forward_solution_sphere(tmp_path, fname_src_small): assert fwd["mri_head_t"]["trans"][0, 3] == -0.05 -@testing.requires_testing_data def test_make_forward_sphere_exclude(): """Test that points are excluded that are outside BEM sphere inner layer.""" - bem = make_sphere_model(r0=(0.0, 0.0, 0.04), head_radius=0.1) + r0 = (0.0, 0.0, 0.04) + inner_radius = 0.08 + head_radius = inner_radius / 0.9 # relative_radii[0] for make_sphere model + bem = make_sphere_model(r0=r0, head_radius=head_radius) + # construct our source space using a sphere the size of the inner (brain) radius, + # with a 1mm exclude zone to avoid numerical issues src = setup_volume_source_space( - "sample", pos=10.0, - sphere=(0.0, 0.0, 0.04, 0.1), - subjects_dir=subjects_dir, + sphere=r0 + (inner_radius,), + mindist=1, # use >0 to avoid any numerical issues exclude=10, - verbose=True, ) - assert src[0]["nuse"] == 3694 + assert src[0]["nuse"] == 2102 # empirically determined trans = Transform("mri", "head") # identity for simplicity raw = read_raw_fif(fname_raw) raw.pick(raw.ch_names[:1]) - fwd = make_forward_solution(raw.info, trans, src, bem, verbose=True) - assert fwd["nsource"] == 3694 - bem_small = make_sphere_model(r0=(0.0, 0.0, 0.04), head_radius=0.09) - fwd_small = make_forward_solution(raw.info, trans, src, bem_small, verbose=True) - assert fwd_small["nsource"] < 3694 + fwd = make_forward_solution(raw.info, trans, src, bem, mindist=0) + assert fwd["nsource"] == src[0]["nuse"] + bem_small = make_sphere_model(r0=r0, head_radius=head_radius - 0.01 / 0.9) + fwd_small = make_forward_solution(raw.info, trans, src, bem_small, mindist=0) + assert fwd_small["nsource"] < src[0]["nuse"] + idx = np.searchsorted(fwd["src"][0]["vertno"], fwd_small["src"][0]["vertno"]) + fwd_data = np.reshape(fwd["sol"]["data"], (len(raw.ch_names), -1, 3)) + fwd_small_data = np.reshape(fwd_small["sol"]["data"], (len(raw.ch_names), -1, 3)) + assert_allclose(fwd_data[:, idx], fwd_small_data) + # again to avoid numerical issues, make this close + fwd_small_2 = make_forward_solution(raw.info, trans, src, bem, mindist=9.999) + assert fwd_small_2["nsource"] == fwd_small["nsource"] + fwd_small_2_data = np.reshape( + fwd_small_2["sol"]["data"], (len(raw.ch_names), -1, 3) + ) + assert_allclose(fwd_small_data, fwd_small_2_data) @pytest.mark.slowtest @@ -879,9 +896,39 @@ def test_sensors_inside_bem(): make_forward_solution(info, trans, fname_src, sphere_noshell) # okay with pytest.raises(RuntimeError, match=".* 42 MEG.*outermost sphere sh.*"): make_forward_solution(info, trans, fname_src, sphere) - sphere = make_sphere_model((0.0, 0.0, 2.0), 1.01) # weird, but okay - make_forward_solution(info, trans, fname_src, sphere) + sphere = make_sphere_model((0.0, 0.0, 0.0), 1.01) # weird, but okay + with pytest.raises(RuntimeError, match=".* 42 MEG.*the outermost sphere shell.*"): + make_forward_solution(info, trans, fname_src, sphere) for ch in info["chs"]: ch["loc"][:3] *= 0.1 with pytest.raises(RuntimeError, match=".* 42 MEG.*the inner skull.*"): make_forward_solution(info, trans, fname_src, fname_bem_meg) + + +def test_make_forward_iterative(): + """Test that points are excluded that are outside BEM sphere inner layer.""" + r0 = (0.0, 0.0, 0.04) + head_radius = 0.08 + bem = make_sphere_model(r0=r0, head_radius=head_radius) + src = setup_volume_source_space( + pos=10.0, + sphere=r0 + (head_radius,), + exclude=10, + ) + assert 500 < src[0]["nuse"] < 4000 + trans = Transform("mri", "head") + raw = read_raw_fif(fname_raw) + raw.pick(raw.ch_names[:2]) + fwd = make_forward_solution(raw.info, trans, src, bem, mindist=0, verbose=True) + # check against iterative version + fm = _ForwardModeler(raw.info, trans, bem) + midpt = fwd["nsource"] // 2 + assert fwd["coord_frame"] == FIFF.FIFFV_COORD_HEAD + fwd_data = list() + rr = apply_trans(invert_transform(fwd["mri_head_t"]), fwd["source_rr"]) + nn = apply_trans(invert_transform(fwd["mri_head_t"]), fwd["source_nn"], move=False) + for sl in (slice(None, midpt), slice(midpt, None)): + ss = setup_volume_source_space(pos=dict(rr=rr[sl], nn=nn[sl])) + fwd_data.append(fm.compute(ss)["sol"]["data"]) + fwd_data = np.concatenate(fwd_data, axis=1) + assert_allclose(fwd_data, fwd["sol"]["data"]) diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 7e2e27d1961..4441502e2d3 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -472,7 +472,7 @@ def __repr__(self): # noqa: D105 @property def _subject(self): - return self[0].get("subject_his_id", None) + return self[0].get("subject_his_id", None) if len(self) else None def __add__(self, other): """Combine source spaces.""" @@ -2587,7 +2587,7 @@ def _filter_source_spaces( if check_inside.inner_r is not None: # ... and those that are at least inner_sphere - limit away mask = ( - np.linalg.norm(check_r1s - check_inside.cm, axis=-1) + np.linalg.norm(check_r1s - check_inside.center, axis=-1) >= check_inside.inner_r - limit / 1000.0 ) idx = idx[mask] diff --git a/mne/surface.py b/mne/surface.py index 23367974e86..34977954cbc 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -760,14 +760,14 @@ def __init__(self, surf, *, mode="old", verbose=None): def _init_old(self): self.inner_r = None - self.cm = self.surf["rr"].mean(0) + self.center = self.surf["rr"].mean(0) # We could use Delaunay or ConvexHull here, Delaunay is slightly slower # to construct but faster to evaluate # See https://stackoverflow.com/questions/16750618/whats-an-efficient-way-to-find-if-a-point-lies-in-the-convex-hull-of-a-point-cl # noqa self.del_tri = Delaunay(self.surf["rr"]) - if self.del_tri.find_simplex(self.cm) >= 0: + if self.del_tri.find_simplex(self.center) >= 0: # Immediately cull some points from the checks - dists = np.linalg.norm(self.surf["rr"] - self.cm, axis=-1) + dists = np.linalg.norm(self.surf["rr"] - self.center, axis=-1) self.inner_r = dists.min() self.outer_r = dists.max() @@ -815,7 +815,7 @@ def _call_old(self, rr, n_jobs): # Limit to indices that can plausibly be outside the surf # but are not definitely outside it if self.inner_r is not None: - dists = np.linalg.norm(rr - self.cm, axis=-1) + dists = np.linalg.norm(rr - self.center, axis=-1) in_mask = dists < self.inner_r n = (in_mask).sum() n_pad = str(n).rjust(prec) @@ -864,12 +864,19 @@ def _call_old(self, rr, n_jobs): class _CheckInsideSphere: - def __init__(self, sphere): + def __init__(self, sphere, *, check="inner"): from .bem import ConductorModel assert isinstance(sphere, ConductorModel) and sphere["is_sphere"] - self.cm = sphere["r0"] - self.inner_r = sphere.radius # float or None + self.center = sphere["r0"] + assert isinstance(check, str) and check in ("inner", "outer"), check + self.check = check + # for a sphere, our closest point and farthest are the same + if len(sphere["layers"]): + self.inner_r = sphere["layers"][0]["rad"] + self.outer_r = sphere["layers"][-1]["rad"] + else: + self.inner_r = self.outer_r = None # No need for verbose dec here because no MNE code is called that would log def __call__(self, rr, *, n_jobs=None, verbose=None): @@ -878,7 +885,11 @@ def __call__(self, rr, *, n_jobs=None, verbose=None): if self.inner_r is None: return np.ones(rr.shape[0], bool) else: - return np.linalg.norm(rr - self.cm, axis=-1) <= self.inner_r + return np.linalg.norm(rr - self.center, axis=-1) <= self._check_r + + @property + def _check_r(self): + return self.inner_r if self.check == "inner" else self.outer_r def query(self, rr): """Return the distance to the sphere surface for each point.""" @@ -888,7 +899,7 @@ def query(self, rr): if self.inner_r is None: dists = np.full(rr.shape[0], np.inf) else: - dists = np.abs(np.linalg.norm(rr - self.cm, axis=-1) - self.inner_r) + dists = np.abs(np.linalg.norm(rr - self.center, axis=-1) - self._check_r) return dists, idx diff --git a/mne/tests/test_dipole.py b/mne/tests/test_dipole.py index f230eaa4256..c52b2283719 100644 --- a/mne/tests/test_dipole.py +++ b/mne/tests/test_dipole.py @@ -389,10 +389,19 @@ def test_accuracy(): ) evoked.pick("meg") evoked.pick([c for c in evoked.ch_names[::4]]) - for rad, perc_90 in zip((0.09, None), (0.002, 0.004)): + for rad, perc_90 in zip((0.095, None), (0.002, 0.004)): bem = make_sphere_model( - "auto", rad, evoked.info, relative_radii=(0.999, 0.998, 0.997, 0.995) + (0.0, 0.0, 0.04), + rad, + evoked.info, + relative_radii=(0.999, 0.998, 0.997, 0.995), ) + if rad is not None: + # These should end up being sorted, and normed by largest relative radius + assert_allclose(bem["layers"][0]["rad"], 0.995 / 0.999 * rad) + assert_allclose(bem["layers"][-1]["rad"], rad) + else: + assert bem["layers"] == [] src = read_source_spaces(fname_src) fwd = make_forward_solution(evoked.info, None, src, bem) From 9f9263405f4a9c23050b9cad8e3f766451c49a2f Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 4 Sep 2025 11:52:25 -0400 Subject: [PATCH 4/8] TST: Full [circle full] From 03ea5560914c479edcc361bcbc4e7620bf7aeb02 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 4 Sep 2025 13:29:21 -0400 Subject: [PATCH 5/8] FIX: Fixes --- examples/simulation/simulate_raw_data.py | 2 -- .../simulated_raw_data_using_subject_anatomy.py | 2 -- mne/chpi.py | 6 +++--- mne/conftest.py | 15 ++++++++++----- mne/minimum_norm/tests/test_inverse.py | 13 +++++++------ mne/simulation/raw.py | 5 +++-- mne/simulation/tests/test_raw.py | 5 +---- mne/source_space/_source_space.py | 4 +++- mne/tests/test_morph.py | 3 ++- tutorials/forward/30_forward.py | 7 ++++--- 10 files changed, 33 insertions(+), 29 deletions(-) diff --git a/examples/simulation/simulate_raw_data.py b/examples/simulation/simulate_raw_data.py index e413a8deb75..0fbefca4480 100644 --- a/examples/simulation/simulate_raw_data.py +++ b/examples/simulation/simulate_raw_data.py @@ -31,8 +31,6 @@ simulate_sparse_stc, ) -print(__doc__) - data_path = sample.data_path() meg_path = data_path / "MEG" / "sample" raw_fname = meg_path / "sample_audvis_raw.fif" diff --git a/examples/simulation/simulated_raw_data_using_subject_anatomy.py b/examples/simulation/simulated_raw_data_using_subject_anatomy.py index 08f30497998..ce6803e2ebe 100644 --- a/examples/simulation/simulated_raw_data_using_subject_anatomy.py +++ b/examples/simulation/simulated_raw_data_using_subject_anatomy.py @@ -26,8 +26,6 @@ import mne from mne.datasets import sample -print(__doc__) - # %% # In this example, raw data will be simulated for the sample subject, so its # information needs to be loaded. This step will download the data if it not diff --git a/mne/chpi.py b/mne/chpi.py index 9aa3894fbbd..1cc168f3e20 100644 --- a/mne/chpi.py +++ b/mne/chpi.py @@ -38,6 +38,7 @@ pick_types, ) from ._fiff.proj import Projection, setup_proj +from .bem import ConductorModel from .channels.channels import _get_meg_system from .cov import compute_whitener, make_ad_hoc_cov from .dipole import _make_guesses @@ -1343,9 +1344,8 @@ def compute_chpi_locs( # Make some location guesses (1 cm grid) R = np.linalg.norm(meg_coils[0], axis=1).min() - guesses = _make_guesses( - dict(R=R, r0=np.zeros(3)), 0.01, 0.0, 0.005, verbose=safe_false - )[0]["rr"] + sphere = ConductorModel(layers=[dict(rad=R)], r0=np.zeros(3), is_sphere=True) + guesses = _make_guesses(sphere, 0.01, 0.0, 0.005, verbose=safe_false)[0]["rr"] logger.info( f"Computing {len(guesses)} HPI location guesses " f"(1 cm grid in a {R * 100:.1f} cm sphere)" diff --git a/mne/conftest.py b/mne/conftest.py index 5555affb3a5..334f1cbe801 100644 --- a/mne/conftest.py +++ b/mne/conftest.py @@ -729,14 +729,16 @@ def _evoked_cov_sphere(_evoked): evoked.pick(evoked.ch_names[::4]) assert len(evoked.ch_names) == 77 cov = mne.read_cov(fname_cov) - sphere = mne.make_sphere_model("auto", "auto", evoked.info) + sphere = mne.make_sphere_model( + (0.0, 0.0, 0.04), 0.1, relative_radii=(0.995, 0.997, 0.998, 1.0) + ) return evoked, cov, sphere @pytest.fixture(scope="session") def _fwd_surf(_evoked_cov_sphere): """Compute a forward for a surface source space.""" - evoked, cov, sphere = _evoked_cov_sphere + evoked, _, sphere = _evoked_cov_sphere src_surf = mne.read_source_spaces(fname_src) return mne.make_forward_solution( evoked.info, fname_trans, src_surf, sphere, mindist=5.0 @@ -747,7 +749,7 @@ def _fwd_surf(_evoked_cov_sphere): def _fwd_subvolume(_evoked_cov_sphere): """Compute a forward for a surface source space.""" pytest.importorskip("nibabel") - evoked, cov, sphere = _evoked_cov_sphere + evoked, _, sphere = _evoked_cov_sphere volume_labels = ["Left-Cerebellum-Cortex", "right-Cerebellum-Cortex"] with pytest.raises(ValueError, match=r"Did you mean one of \['Right-Cere"): mne.setup_volume_source_space( @@ -761,9 +763,12 @@ def _fwd_subvolume(_evoked_cov_sphere): subjects_dir=subjects_dir, add_interpolator=False, ) - return mne.make_forward_solution( - evoked.info, fname_trans, src_vol, sphere, mindist=5.0 + fwd = mne.make_forward_solution( + evoked.info, fname_trans, src_vol, sphere, mindist=1.0 ) + nsrc = sum(s["nuse"] for s in src_vol) + assert fwd["nsource"] == nsrc + return fwd @pytest.fixture diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 5b5c941a9ac..88578b20023 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -543,7 +543,7 @@ def test_apply_inverse_sphere(evoked, tmp_path): with evoked.info._unlock(): evoked.info["projs"] = [] cov = make_ad_hoc_cov(evoked.info) - sphere = make_sphere_model("auto", "auto", evoked.info) + sphere = make_sphere_model((0.0, 0.0, 0.04), None) fwd = read_forward_solution(fname_fwd) vertices = [fwd["src"][0]["vertno"][::5], fwd["src"][1]["vertno"][::5]] stc = SourceEstimate( @@ -1577,11 +1577,11 @@ def test_inverse_mixed_loose(mixed_fwd_cov_evoked): evoked_sim = EvokedArray(data, evoked.info) del data # dipole - sphere = mne.make_sphere_model("auto", "auto", evoked.info) + sphere = mne.make_sphere_model((0.0, 0.0, 0.05), 0.1) dip, _ = mne.fit_dipole(evoked_sim, cov, sphere) assert_allclose(dip.pos, want_pos, atol=1e-2) # 1 cm ang = np.rad2deg(np.arccos(np.sum(dip.ori * want_ori, axis=1))) - assert_array_less(ang, 65) # not great + assert_array_less(ang, 70) # not great # MNE stc = apply_inverse(evoked_sim, inv_fixed, pick_ori="vector") stc, nn = stc.project("pca", fwd["src"]) @@ -1594,9 +1594,10 @@ def test_inverse_mixed_loose(mixed_fwd_cov_evoked): ) got_ori = nn[idx] got_pos = fwd["source_rr"][idx] - assert_allclose(got_pos, want_pos, atol=1.1e-2) # 1.1 cm + assert_allclose(got_pos, want_pos, atol=3e-2) # 3 cm ang = np.rad2deg(np.arccos(np.sum(got_ori * want_ori, axis=1))) - assert_array_less(ang, 40) # better than ECD + sphere + ang = np.minimum(ang, 180 - ang) # don't care about direction + assert_array_less(ang, 50) # better than ECD + sphere # MxNE stc = mne.inverse_sparse.mixed_norm( evoked, @@ -1613,7 +1614,7 @@ def test_inverse_mixed_loose(mixed_fwd_cov_evoked): assert len(stc.data) == 2 pos = np.concatenate([fwd["src"][ii]["rr"][v] for ii, v in enumerate(stc.vertices)]) assert pos.shape == (2, 3) - assert_allclose(got_pos, want_pos, atol=1.1e-2) + assert_allclose(got_pos, want_pos, atol=3e-2) @testing.requires_testing_data diff --git a/mne/simulation/raw.py b/mne/simulation/raw.py index 298cd9bf185..94f186992e1 100644 --- a/mne/simulation/raw.py +++ b/mne/simulation/raw.py @@ -502,9 +502,10 @@ def _add_exg(raw, kind, head_pos, interp, n_jobs, random_state): meg_picks = pick_types(info, meg=True, eeg=False, exclude=()) meeg_picks = pick_types(info, meg=True, eeg=True, exclude=()) R, r0 = fit_sphere_to_headshape(info, units="m", verbose=False)[:2] + head_radius = R if kind == "blink" else None bem = make_sphere_model( r0, - head_radius=R, + head_radius=head_radius, relative_radii=(0.97, 0.98, 0.99, 1.0), sigmas=(0.33, 1.0, 0.004, 0.33), verbose=False, @@ -849,7 +850,7 @@ def _iter_forward_solutions( # Compute forward if forward is None: if not bem["is_sphere"]: - outside = ~_CheckInside(bem_surf)(coil_rr, n_jobs, verbose=False) + outside = ~_CheckInside(bem_surf)(coil_rr, n_jobs=n_jobs, verbose=False) elif bem.radius is not None: d = coil_rr - bem["r0"] outside = np.sqrt(np.sum(d * d, axis=1)) > bem.radius diff --git a/mne/simulation/tests/test_raw.py b/mne/simulation/tests/test_raw.py index 08e05ed6eca..69194681467 100644 --- a/mne/simulation/tests/test_raw.py +++ b/mne/simulation/tests/test_raw.py @@ -223,10 +223,7 @@ def raw_data(): src = read_source_spaces(src_fname) trans = read_trans(trans_fname) # Use fixed values from old sphere fit to reduce lines changed with fixed algorithm - sphere = make_sphere_model( - [-0.00413508, 0.01598787, 0.05175598], - 0.09100286249131773, - ) + sphere = make_sphere_model([0.0, 0.0, 0.03], 0.1) stc = _make_stc(raw, src) return raw, src, stc, trans, sphere diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 4441502e2d3..42458151af9 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -2142,7 +2142,9 @@ def _make_volume_source_space( kwargs = dict(limit=mindist, mri_head_t=None, src=[sp]) assert sp["coord_frame"] == FIFF.FIFFV_COORD_MRI if "rr" in surf: - assert surf["coord_frame"] == FIFF.FIFFV_COORD_MRI + # If it's a surface passed as a dict, it might not have a coord_frame, + # assume it's been read by read_surface and is in mri coords + assert surf.get("coord_frame", FIFF.FIFFV_COORD_MRI) == FIFF.FIFFV_COORD_MRI else: assert surf["is_sphere"] _filter_source_spaces(surf, n_jobs=n_jobs, **kwargs) diff --git a/mne/tests/test_morph.py b/mne/tests/test_morph.py index 7cb7d0cb9d9..9296e13ce39 100644 --- a/mne/tests/test_morph.py +++ b/mne/tests/test_morph.py @@ -920,8 +920,9 @@ def test_volume_labels_morph(tmp_path, sl, n_real, n_mri, n_orig): assert len(src) == n_use assert src.kind == "volume" n_src = sum(s["nuse"] for s in src) - sphere = make_sphere_model("auto", "auto", evoked.info) + sphere = make_sphere_model((0.0, 0.0, 0.04), 0.1) fwd = make_forward_solution(evoked.info, fname_trans, src, sphere) + assert fwd["nsource"] == n_src assert fwd["sol"]["data"].shape == (n_ch, n_src * 3) inv = make_inverse_operator( evoked.info, fwd, make_ad_hoc_cov(evoked.info), loose=1.0 diff --git a/tutorials/forward/30_forward.py b/tutorials/forward/30_forward.py index 72731982962..ef39786cdf5 100644 --- a/tutorials/forward/30_forward.py +++ b/tutorials/forward/30_forward.py @@ -25,6 +25,7 @@ # the raw file containing the channel location + types sample_dir = data_path / "MEG" / "sample" raw_fname = sample_dir / "sample_audvis_raw.fif" + # The paths to Freesurfer reconstructions subjects_dir = data_path / "subjects" subject = "sample" @@ -35,9 +36,9 @@ # # To compute a forward operator we need: # -# - a ``-trans.fif`` file that contains the coregistration info. -# - a source space -# - the :term:`BEM` surfaces +# - a ``-trans.fif`` file that contains the coregistration info. +# - a source space +# - the :term:`BEM` surfaces # %% # Compute and visualize BEM surfaces From a32ea0999c98fb4cb1005cf31263673f8a8591b4 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Fri, 5 Sep 2025 11:28:00 -0400 Subject: [PATCH 6/8] FIX: Tests --- mne/forward/tests/test_make_forward.py | 6 +++++- mne/inverse_sparse/tests/test_mxne_inverse.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index a82b4cb150b..34494c57fde 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -552,7 +552,11 @@ def test_make_forward_solution_sphere(tmp_path, fname_src_small): ] ) fwd = read_forward_solution(out_name) - sphere = make_sphere_model(verbose=True) + sphere = make_sphere_model( + head_radius=0.1, + relative_radii=(0.95, 0.97, 0.98, 1), + verbose=True, + ) src = read_source_spaces(fname_src_small) fwd_py = make_forward_solution( fname_raw, fname_trans, src, sphere, meg=True, eeg=True, verbose=True diff --git a/mne/inverse_sparse/tests/test_mxne_inverse.py b/mne/inverse_sparse/tests/test_mxne_inverse.py index cff8adc3225..4a3a15df0f4 100644 --- a/mne/inverse_sparse/tests/test_mxne_inverse.py +++ b/mne/inverse_sparse/tests/test_mxne_inverse.py @@ -306,7 +306,7 @@ def test_mxne_vol_sphere(): evoked_l21.crop(tmin=0.081, tmax=0.1) info = evoked.info - sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.080) + sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.09) src = mne.setup_volume_source_space( subject=None, pos=15.0, From 1610f6137cdc17445df7569066106fcc5f29afd7 Mon Sep 17 00:00:00 2001 From: Daniel McCloy Date: Fri, 5 Sep 2025 13:25:35 -0500 Subject: [PATCH 7/8] fix for pip-pre --- mne/_ola.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/_ola.py b/mne/_ola.py index e43e7cd3d31..939aff39646 100644 --- a/mne/_ola.py +++ b/mne/_ola.py @@ -309,7 +309,7 @@ def __init__( # Create our window boundaries window_name = window if isinstance(window, str) else "custom" self._window = get_window( - window, self._n_samples, fftbins=(self._n_samples - 1) % 2 + window, self._n_samples, fftbins=bool((self._n_samples - 1) % 2) ) self._window /= _check_cola( self._window, self._n_samples, self._step, window_name, tol=tol From 1a9899e7e551be0cfe1d16538cb6521ceee1f1d7 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 8 Sep 2025 10:30:59 -0400 Subject: [PATCH 8/8] FIX: DRY --- mne/forward/_make_forward.py | 9 ++------- mne/forward/forward.py | 27 ++++++--------------------- mne/minimum_norm/inverse.py | 13 ++----------- mne/simulation/tests/test_raw.py | 4 +--- mne/source_space/_source_space.py | 11 +++++++++++ mne/tests/test_dipole.py | 4 +--- mne/tests/test_source_estimate.py | 5 ++--- 7 files changed, 25 insertions(+), 48 deletions(-) diff --git a/mne/forward/_make_forward.py b/mne/forward/_make_forward.py index 7342574b704..02eb6cc1c6a 100644 --- a/mne/forward/_make_forward.py +++ b/mne/forward/_make_forward.py @@ -36,7 +36,6 @@ _print_coord_trans, apply_trans, invert_transform, - transform_surface_to, ) from ..utils import _check_fname, _pl, _validate_type, logger, verbose, warn from ._compute_forward import ( @@ -522,9 +521,7 @@ def _prepare_for_forward( # Transform the source spaces into the appropriate coordinates # (will either be HEAD or MRI) - for s in src: - transform_surface_to(s, "head", mri_head_t) - del s + src._transform_to("head", mri_head_t) if len(src): logger.info( f"Source spaces are now in {_coord_frame_name(src[0]['coord_frame'])} " @@ -990,9 +987,7 @@ def __init__( def compute(self, src): src = _ensure_src(src).copy() - for s in src: - transform_surface_to(s, "head", self.mri_head_t) - + src._transform_to("head", self.mri_head_t) kwargs = dict(limit=self.mindist, mri_head_t=self.mri_head_t, src=src) _filter_source_spaces(self.check_inside, n_jobs=self.n_jobs, **kwargs) rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) diff --git a/mne/forward/forward.py b/mne/forward/forward.py index f014b69a64e..e8a7d62a3ce 100644 --- a/mne/forward/forward.py +++ b/mne/forward/forward.py @@ -62,7 +62,7 @@ find_source_space_hemi, ) from ..surface import _normal_orth -from ..transforms import invert_transform, transform_surface_to, write_trans +from ..transforms import invert_transform, write_trans from ..utils import ( _check_compensation_grade, _check_fname, @@ -680,16 +680,8 @@ def read_forward_solution(fname, include=(), exclude=(), *, ordered=True, verbos # Transform each source space to the HEAD or MRI coordinate frame, # depending on the coordinate frame of the forward solution - # NOTE: the function transform_surface_to will also work on discrete and - # volume sources - nuse = 0 - for s in src: - try: - s = transform_surface_to(s, fwd["coord_frame"], mri_head_t) - except Exception as inst: - raise ValueError(f"Could not transform source space ({inst})") - - nuse += s["nuse"] + src._transform_to(fwd["coord_frame"], mri_head_t) + nuse = sum(s["nuse"] for s in src) # Make sure the number of sources match after transformation if nuse != fwd["nsource"]: @@ -952,16 +944,9 @@ def _write_forward_solution(fid, fwd): write_forward_meas_info(fid, fwd["info"]) # invert our original source space transform - src = list() - for s in fwd["src"]: - s = deepcopy(s) - try: - # returns source space to original coordinate frame - # usually MRI - s = transform_surface_to(s, fwd["mri_head_t"]["from"], fwd["mri_head_t"]) - except Exception as inst: - raise ValueError(f"Could not transform source space ({inst})") - src.append(s) + src = fwd["src"].copy() + # returns source space to original coordinate frame, usually MRI + src._transform_to(fwd["mri_head_t"]["from"], fwd["mri_head_t"]) # # Write the source spaces (again) diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index 7c789503ac1..da770744881 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -62,7 +62,7 @@ ) from ..surface import _normal_orth from ..time_frequency.tfr import _check_tfr_complex -from ..transforms import _ensure_trans, transform_surface_to +from ..transforms import _ensure_trans from ..utils import ( _check_compensation_grade, _check_depth, @@ -381,16 +381,7 @@ def read_inverse_operator(fname, *, verbose=None): inv["reginv"] = [] inv["noisenorm"] = [] # These are the noise-normalization factors # - nuse = 0 - for k in range(len(inv["src"])): - try: - inv["src"][k] = transform_surface_to( - inv["src"][k], inv["coord_frame"], mri_head_t - ) - except Exception as inst: - raise Exception(f"Could not transform source space ({inst})") - - nuse += inv["src"][k]["nuse"] + inv["src"]._transform_to(inv["coord_frame"], mri_head_t) logger.info( " Source spaces transformed to the inverse solution coordinate frame" diff --git a/mne/simulation/tests/test_raw.py b/mne/simulation/tests/test_raw.py index 69194681467..4b1d514e8d7 100644 --- a/mne/simulation/tests/test_raw.py +++ b/mne/simulation/tests/test_raw.py @@ -28,7 +28,6 @@ read_trans, setup_source_space, setup_volume_source_space, - transform_surface_to, ) from mne._fiff.constants import FIFF from mne.bem import _surfaces_to_bem @@ -388,8 +387,7 @@ def test_simulate_raw_bem(raw_data): med_corr = np.median(np.diag(corr[:n_ch, -n_ch:])) assert med_corr > 0.65 # do some round-trip localization - for s in src: - transform_surface_to(s, "head", trans) + src._transform_to("head", trans) locs = np.concatenate([s["rr"][s["vertno"]] for s in src]) tmax = (len(locs) - 1) / raw.info["sfreq"] cov = make_ad_hoc_cov(raw.info) diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 42458151af9..948e7afe8e6 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -68,6 +68,7 @@ apply_trans, combine_transforms, invert_transform, + transform_surface_to, ) from ..utils import ( _check_fname, @@ -787,6 +788,16 @@ def export_volume( # write image to file nib.save(img, fname) + def _transform_to(self, cf, mri_head_t): + # NOTE: the function transform_surface_to will also work on discrete and + # volume sources + try: + for s in self: + transform_surface_to(s, cf, mri_head_t, copy=False) + except Exception as inst: + raise RuntimeError(f"Could not transform source space ({inst})") + return self + def _add_patch_info(s): """Patch information in a source space. diff --git a/mne/tests/test_dipole.py b/mne/tests/test_dipole.py index c52b2283719..c6235171fa3 100644 --- a/mne/tests/test_dipole.py +++ b/mne/tests/test_dipole.py @@ -30,7 +30,6 @@ read_evokeds, read_forward_solution, read_source_spaces, - transform_surface_to, write_evokeds, ) from mne._fiff.constants import FIFF @@ -219,8 +218,7 @@ def test_dipole_fitting(tmp_path): ) # Compare to original points - transform_surface_to(fwd["src"][0], "head", fwd["mri_head_t"]) - transform_surface_to(fwd["src"][1], "head", fwd["mri_head_t"]) + fwd["src"]._transform_to("head", fwd["mri_head_t"]) assert fwd["src"][0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD src_rr = np.concatenate([s["rr"][v] for s, v in zip(fwd["src"], vertices)], axis=0) src_nn = np.concatenate([s["nn"][v] for s, v in zip(fwd["src"], vertices)], axis=0) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index b64ba792a6c..41e33c52274 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -75,7 +75,7 @@ from mne.morph_map import _make_morph_map_hemi from mne.source_estimate import _get_vol_mask, _make_stc, grade_to_tris from mne.source_space._source_space import _get_src_nn -from mne.transforms import apply_trans, invert_transform, transform_surface_to +from mne.transforms import apply_trans, invert_transform from mne.utils import ( _record_warnings, catch_logging, @@ -1709,8 +1709,7 @@ def test_stc_near_sensors(tmp_path): assert np.isclose(dists, 0.0, atol=1e-6).any(0).all() src = read_source_spaces(fname_src) # uses "white" but should be okay - for s in src: - transform_surface_to(s, "head", trans, copy=False) + src._transform_to("head", trans) assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD stc_src = stc_near_sensors(evoked, src=src, **kwargs) assert len(stc_src.data) == 7928