diff --git a/doc/api/inverse.rst b/doc/api/inverse.rst index 754244c17fe..9ff9e65655f 100644 --- a/doc/api/inverse.rst +++ b/doc/api/inverse.rst @@ -85,6 +85,7 @@ Inverse Solutions Dipole DipoleFixed fit_dipole + gui.dipolefit :py:mod:`mne.dipole`: diff --git a/doc/changes/dev/13074.newfeature.rst b/doc/changes/dev/13074.newfeature.rst new file mode 100644 index 00000000000..35415ea90d7 --- /dev/null +++ b/doc/changes/dev/13074.newfeature.rst @@ -0,0 +1 @@ +Add a GUI for interactive guided dipole fitting (:func:`mne.gui.dipolefit`), by `Marijn van Vliet`_ diff --git a/doc/conf.py b/doc/conf.py index f4eb8eda8f6..487e86ab884 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -442,6 +442,7 @@ "_Renderer", "n_triangles", "CoregistrationUI", + "DipoleFitUI", "mne_qt_browser.figure.MNEQtBrowser", # pooch, since its website is unreliable and users will rarely need the links "pooch.Unzip", diff --git a/doc/sphinxext/mne_doc_utils.py b/doc/sphinxext/mne_doc_utils.py index 1c419aaa18c..3afc1e30cd8 100644 --- a/doc/sphinxext/mne_doc_utils.py +++ b/doc/sphinxext/mne_doc_utils.py @@ -144,6 +144,12 @@ def reset_modules(gallery_conf, fname, when): """Do the reset.""" import matplotlib.pyplot as plt + # Examples that set ``# sphinx_gallery_preserve_gui = True`` keep a single GUI open + # across all of their code blocks, and the scraper (not the example's globals, which + # sphinx-gallery has already dropped by the time we get here with when="after") + # holds the last reference to it. Close them before the leak checks below. + gui_scraper.close_preserved() + mne.viz.set_3d_backend("pyvistaqt") pyvista.OFF_SCREEN = False pyvista.BUILDING_GALLERY = True diff --git a/mne/commands/mne_dipolefit.py b/mne/commands/mne_dipolefit.py new file mode 100644 index 00000000000..e6b8b6224e4 --- /dev/null +++ b/mne/commands/mne_dipolefit.py @@ -0,0 +1,166 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +"""Open the dipole fitting GUI on the given evoked file ("-ave.fif"). + +Examples +-------- +.. code-block:: console + + $ mne dipolefit + +""" + +import os.path as op +import sys + +import mne + + +def run(): + """Run command.""" + from mne.commands.utils import _add_verbose_flag, get_optparser + + parser = get_optparser(__file__, usage="mne dipolefit EVOKED_FILE") + + parser.add_option( + "--condition", + default=0, + help="The condition to use.", + ) + parser.add_option( + "--baseline", + default=None, + metavar="BASELINE", # docutils doesn't like commas in rendered metavar + help=( + "The time period to use as baseline, written as two numbers (in seconds, " + "relative to the stimulus onset) separated by a comma. " + "For example: --baseline=-0.2,0.1" + ), + ) + parser.add_option( + "-c", + "--cov", + default=None, + metavar="COV_FILE", + help='The noise covariance ("-cov.fif") to use.', + ) + parser.add_option( + "-b", + "--bem", + default=None, + metavar="BEM_FILE", + help=( + 'The BEM model ("-bem-sol.fif") to use. When omitted, a basic sphere ' + "model will be used." + ), + ) + parser.add_option( + "-t", + "--initial-time", + default=None, + type=float, + metavar="TIME", + help="The initial time to show", + ) + parser.add_option( + "--trans", + default=None, + metavar="TRANS_FILE", + help='Head<->MRI transform FIF file ("-trans.fif")', + ) + parser.add_option( + "--stc", + default=None, + metavar="STC_FILE", + help="An optional distributed source estimate to show during dipole fitting.", + ) + parser.add_option( + "-s", "--subject", dest="subject", default=None, help="Subject name" + ) + parser.add_option( + "-d", + "--subjects-dir", + default=None, + help="Subjects directory", + ) + parser.add_option( + "--hide-density", + action="store_true", + default=False, + help="Prevent showing the magnetic field density as blobs of color.", + ) + parser.add_option( + "--channel-type", + default=None, + help=( + 'Restrict channel types to either "meg" or "eeg". By default both are used ' + "if present." + ), + ) + parser.add_option( + "-j", "--n-jobs", default=-1, type=int, help="Number of CPUs to use." + ) + _add_verbose_flag(parser) + + options, args = parser.parse_args() + if len(args) != 1: + parser.print_help() + sys.exit(1) + + # expanduser allows ~ for paths + subjects_dir = options.subjects_dir + if subjects_dir is not None: + subjects_dir = op.expanduser(subjects_dir) + bem = options.bem + if bem is not None: + bem = op.expanduser(bem) + trans = options.trans + if trans is not None: + trans = op.expanduser(trans) + stc = options.stc + if stc is not None: + stc = op.expanduser(stc) + + # Condition can be specified as integer index or string comment. + try: + condition = int(options.condition) + except ValueError: + condition = options.condition + evoked = mne.read_evokeds(args[0], condition=condition) + + # Parse the baseline time period + baseline = None + if options.baseline: + try: + baseline = [float(x) for x in options.baseline.split(",")] + if len(baseline) != 2: + raise ValueError() + except ValueError: + raise ValueError( + "The 'baseline' parameter should be written as two numbers (in seconds," + " relative to the stimulus onset) separated by a comma. " + "For example: --baseline=-0.2,0.1" + ) + + mne.gui.dipolefit( + evoked=evoked, + baseline=baseline, + cov=options.cov, + bem=bem, + subject=options.subject, + subjects_dir=subjects_dir, + stc=stc, + ch_type=options.channel_type, + initial_time=options.initial_time, + trans=trans, + n_jobs=options.n_jobs, + show_density=not options.hide_density, + show=True, + block=True, + verbose=options.verbose, + ) + + +mne.utils.run_command_if_main() diff --git a/mne/commands/tests/test_commands.py b/mne/commands/tests/test_commands.py index 41347cec21b..27aa3878011 100644 --- a/mne/commands/tests/test_commands.py +++ b/mne/commands/tests/test_commands.py @@ -31,6 +31,7 @@ mne_compute_proj_ecg, mne_compute_proj_eog, mne_coreg, + mne_dipolefit, mne_flash_bem, mne_kit2fiff, mne_make_scalp_surfaces, @@ -60,9 +61,9 @@ raw_fname = op.join(base_dir, "test_raw.fif") testing_path = testing.data_path(download=False) -subjects_dir = op.join(testing_path, "subjects") -bem_model_fname = op.join( - testing_path, "subjects", "sample", "bem", "sample-320-320-320-bem.fif" +subjects_dir = testing_path / "subjects" +bem_model_fname = ( + testing_path / "subjects" / "sample" / "bem" / "sample-320-320-320-bem.fif" ) @@ -619,3 +620,53 @@ def test_anonymize(tmp_path): info = read_info(out_fname) assert op.exists(out_fname) assert info["meas_date"] == _stamp_to_dt((946684800, 0)) + + +def test_dipolefit(monkeypatch): + """Test mne dipolefit.""" + check_usage(mne_dipolefit) + # Don't open the GUI, just check that the arguments are passed along correctly. + kwargs = dict() + monkeypatch.setattr(mne.gui, "dipolefit", lambda **kw: kwargs.update(kw)) + ave_fname = op.join(base_dir, "test-ave.fif") + args = ( + ave_fname, + "--condition=Right Auditory", + "--baseline=-0.2,0", + "--channel-type=meg", + "--initial-time=0.1", + "--hide-density", + "--subject=fake", + "--subjects-dir=~/fake-subjects", + "--bem=~/fake-bem-sol.fif", + "--trans=~/fake-trans.fif", + "--stc=~/fake-stc", + ) + with ArgvSetter(args): + mne_dipolefit.run() + assert kwargs["evoked"].comment == "Right Auditory" + assert kwargs["baseline"] == [-0.2, 0] + assert kwargs["ch_type"] == "meg" + assert kwargs["initial_time"] == 0.1 + assert kwargs["show_density"] is False + assert kwargs["subject"] == "fake" + for key, val in dict( + subjects_dir="~/fake-subjects", + bem="~/fake-bem-sol.fif", + trans="~/fake-trans.fif", + stc="~/fake-stc", + ).items(): + assert kwargs[key] == op.expanduser(val) # "~" gets expanded + + # The condition can also be given as an index (the default being the first one). + with ArgvSetter((ave_fname,)): + mne_dipolefit.run() + assert kwargs["evoked"].comment == "Left Auditory" + assert kwargs["baseline"] is None + assert kwargs["show_density"] is True + assert kwargs["bem"] is kwargs["trans"] is kwargs["stc"] is None + + # The baseline needs to be two comma-separated numbers. + with ArgvSetter((ave_fname, "--baseline=0")): + with pytest.raises(ValueError, match="two numbers"): + mne_dipolefit.run() diff --git a/mne/dipole.py b/mne/dipole.py index 3833491be75..d83fabd9333 100644 --- a/mne/dipole.py +++ b/mne/dipole.py @@ -914,7 +914,7 @@ def _write_dipole_bdip(fname, dip): fid.write(np.array(has_errors, ">i4").tobytes()) # has_errors fid.write(np.zeros(1, ">f4").tobytes()) # noise level for key in _BDIP_ERROR_KEYS: - val = dip.conf[key][ti] if key in dip.conf else 0.0 + val = dip.conf[key][ti] if key in dip.conf else np.array(0.0) assert val.shape == () fid.write(np.array(val, ">f4").tobytes()) fid.write(np.zeros(25, ">f4").tobytes()) diff --git a/mne/gui/__init__.pyi b/mne/gui/__init__.pyi index 086c51a4904..a6dc001ef84 100644 --- a/mne/gui/__init__.pyi +++ b/mne/gui/__init__.pyi @@ -1,2 +1,2 @@ -__all__ = ["_GUIScraper", "coregistration"] -from ._gui import _GUIScraper, coregistration +__all__ = ["_GUIScraper", "coregistration", "dipolefit"] +from ._gui import _GUIScraper, coregistration, dipolefit diff --git a/mne/gui/_dipolefit.py b/mne/gui/_dipolefit.py new file mode 100644 index 00000000000..7e138e258bb --- /dev/null +++ b/mne/gui/_dipolefit.py @@ -0,0 +1,985 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +from copy import deepcopy +from functools import partial +from pathlib import Path + +import numpy as np +import pyvista + +from .._fiff.pick import pick_types +from ..bem import ( + ConductorModel, + _ensure_bem_surfaces, + make_sphere_model, + read_bem_solution, +) +from ..cov import _ensure_cov, make_ad_hoc_cov +from ..dipole import Dipole, fit_dipole +from ..evoked import Evoked +from ..forward import convert_forward_solution, make_field_map +from ..forward._make_forward import _ForwardModeler +from ..minimum_norm import apply_inverse, make_inverse_operator +from ..source_estimate import ( + SourceEstimate, + _BaseSurfaceSourceEstimate, + read_source_estimate, +) +from ..source_space import setup_volume_source_space +from ..surface import _normal_orth +from ..transforms import _get_trans, _get_transforms_to_coord_frame, apply_trans +from ..utils import ( + _auto_weakref, + _check_option, + _validate_type, + fill_doc, + logger, + verbose, +) +from ..viz import EvokedField, create_3d_figure +from ..viz._3d import _plot_head_surface, _plot_sensors_3d +from ..viz.backends._utils import _qt_app_exec +from ..viz.ui_events import ChannelsSelect, TimeChange, link, publish, subscribe +from ..viz.utils import _get_color_list + + +@fill_doc +class DipoleFitUI: + """GUI for interactive dipole fitting, inspired by MEGIN's XFit program. + + Parameters + ---------- + evoked : instance of Evoked | path-like + Evoked data to show fieldmap of and fit dipoles to. + %(baseline_evoked)s + cov : instance of Covariance | "baseline" | None + Noise covariance matrix. If ``None``, an ad-hoc covariance matrix is used with + default values for the diagonal elements (see Notes). If ``"baseline"``, the + diagonal elements is estimated from the baseline period of the evoked data. + bem : instance of ConductorModel | path-like | None + Boundary element model to use in forward calculations, or a path to the BEM + solution file (``"-bem-sol.fif"``) to read it from. If ``None``, a spherical + model is used. + initial_time : float | None + Initial time point to show. If ``None``, the time point of the maximum field + strength is used. + trans : instance of Transform | None + The transformation from head coordinates to MRI coordinates. If ``None``, + the identity matrix is used and everything will be done in head coordinates. + stc : instance of SourceEstimate | None + An optional distributed source estimate to show alongside the fieldmap. The time + samples need to match those of the evoked data. + subject : str | None + The subject name. If ``None``, no MRI data is shown. + %(subjects_dir)s + surf_maps : list | None + The surface mapping information obtained with make_field_map. If ``None``, one + will be generated based on the given data. + %(rank)s + show_density : bool + Whether to show the density of the fieldmap. + ch_type : "meg" | "eeg" | None + Type of channels to use for the dipole fitting. By default (``None``) both MEG + and EEG channels will be used. + show_sensors : bool + Whether to show the sensors in the 3D view. + %(n_jobs)s + show : bool + Show the GUI if True. + block : bool + Whether to halt program execution until the figure is closed. + %(verbose)s + + Attributes + ---------- + dipoles : list of Dipole + All currently enabled dipoles in the model. + """ + + def __init__( + self, + evoked, + *, + baseline=None, + cov=None, + bem=None, + initial_time=None, + trans=None, + stc=None, + subject=None, + subjects_dir=None, + surf_maps=None, + rank="info", + show_density=True, + ch_type=None, + show_sensors=True, + n_jobs=None, + show=True, + block=False, + verbose=None, + ): + _validate_type(evoked, Evoked, "evoked") + if baseline is not None: + evoked = evoked.copy().apply_baseline(baseline) + + if cov is None: + logger.info("Using ad-hoc noise covariance.") + cov = make_ad_hoc_cov(evoked.info) + elif cov == "baseline": + if evoked.baseline is None: + raise ValueError( + 'cov="baseline" requires baseline-corrected data. Set the ' + "baseline parameter or baseline-correct the evoked data first." + ) + logger.info( + f"Estimating noise covariance from baseline ({evoked.baseline[0]:.3f} " + f"to {evoked.baseline[1]:.3f} seconds)." + ) + std = dict() + for typ in set(evoked.get_channel_types(only_data_chs=True)): + baseline = evoked.copy().pick(typ).crop(*evoked.baseline) + std[typ] = baseline.data.std(axis=1).mean() + cov = make_ad_hoc_cov(evoked.info, std) + else: + cov = _ensure_cov(cov) + + _validate_type(bem, ("path-like", ConductorModel, None), "bem") + if bem is None: + bem = make_sphere_model("auto", "auto", evoked.info) + elif not isinstance(bem, ConductorModel): + # a path means a BEM solution file (cf. _make_forward._setup_bem) + bem = read_bem_solution(bem) + bem = _ensure_bem_surfaces(bem, extra_allow=(ConductorModel,)) + + if ch_type is not None: + evoked = evoked.copy().pick(ch_type) + + if surf_maps is None: + surf_maps = make_field_map( + evoked, + trans=trans, + origin=bem["r0"] if bem["is_sphere"] else "auto", + subject=subject, + subjects_dir=subjects_dir, + n_jobs=n_jobs, + verbose=verbose, + ) + + if initial_time is None: + # Set initial time to moment of maximum field power. + data = evoked.copy().pick(surf_maps[0]["ch_names"]).data + initial_time = evoked.times[np.argmax(np.mean(data**2, axis=0))] + + if stc is not None: + _validate_type(stc, ("path-like", _BaseSurfaceSourceEstimate), "stc") + if not isinstance(stc, _BaseSurfaceSourceEstimate): + stc = read_source_estimate(stc) + + if len(stc.times) != len(evoked.times) or not np.allclose( + stc.times, evoked.times + ): + raise ValueError( + "The time samples of the source estimate do not match those of the " + "evoked data." + ) + if trans is None: + raise ValueError( + "`trans` cannot be `None` when showing the fieldlines in " + "combination with a source estimate." + ) + + # Get transforms to convert all the various meshes to MRI space. + head_mri_t = _get_trans(trans, "head", "mri")[0] + to_cf_t = _get_transforms_to_coord_frame( + evoked.info, head_mri_t, coord_frame="mri" + ) + + self.fwd = _ForwardModeler( + info=evoked.info, + trans=trans, + bem=bem, + n_jobs=n_jobs, + verbose=verbose, + ) + + # Initialize all the private attributes. + self._actors = dict() + self._bem = bem + self._ch_type = ch_type + self._cov = cov + self._current_time = initial_time + self._dipoles = dict() + self._evoked = evoked + self._helmet_surf = None + self._surf_maps = surf_maps + self._fig_sensors = None + self._multi_dipole_method = "Multi dipole (MNE)" + self._show_density = show_density + self._stc = stc + self._subjects_dir = subjects_dir + self._subject = subject + self._time_line = None + self._head_mri_t = head_mri_t + self._to_cf_t = to_cf_t + self._rank = rank + self._verbose = verbose + self._n_jobs = n_jobs + + # Configure the GUI. + self._configure_main_display( + show_sensors=show_sensors, show=show + ) # sets self._fig + self._configure_dock() + + # must be done last + if show: + self._renderer.show() + if block and self._renderer._kind != "notebook": + _qt_app_exec(self._renderer.figure.store["app"]) + + @property + def _renderer(self): + return self._fig._renderer + + @property + def dipoles(self): + """A list of all the fitted dipoles that are enabled in the GUI.""" + return [d["dip"] for d in self._dipoles.values() if d["active"]] + + def _configure_main_display(self, show_sensors=True, show=True): + """Configure main 3D display of the GUI.""" + fig_into = create_3d_figure((1080, 720), bgcolor="white", show=show) + + self._stc_brain = None + if self._stc is not None: + kwargs = dict( + subject=self._subject, + subjects_dir=self._subjects_dir, + hemi="both", + time_viewer=False, + initial_time=self._current_time, + brain_kwargs=dict(units="m"), + figure=fig_into, + ) + if isinstance(self._stc, SourceEstimate): + kwargs["surface"] = "white" + self._stc_brain = self._stc.plot(**kwargs) + self._actors["brain"] = self._stc_brain._actors["data"] + fig_into = self._stc_brain # plot into the brain instead + + fig_ef = EvokedField( + self._evoked, + self._surf_maps, + time=self._current_time, + interpolation="linear", + alpha=0, + show_density=self._show_density, + foreground="black", + background="white", + fig=fig_into, # can be Figure3D or Brain instance + ) + del fig_into + fig_ef.separate_canvas = False # needed to plot the timeline later + fig_ef.set_contour_line_width(2) + if self._stc is not None: + link(self._stc_brain, fig_ef) + + for surf_map in fig_ef._surf_maps: + if surf_map["map_kind"] == "meg": + helmet_mesh = surf_map["mesh"] + helmet_mesh._actor.prop.culling = "back" + self._actors["helmet"] = helmet_mesh._actor + # needed later to draw the big arrows on the helmet + self._helmet_surf = surf_map["surf"] + # For MEG fieldlines, we want to occlude the ones not facing us, + # otherwise it's hard to interpret them. Since the "contours" object + # does not support backface culling, we create an opaque mesh to put in + # front of the contour lines with frontface culling. + occl_surf = deepcopy(surf_map["surf"]) + occl_surf["rr"] -= 1e-3 * occl_surf["nn"] + occl_act, _ = fig_ef._renderer.surface(occl_surf, color="white") + occl_act.prop.culling = "front" + occl_act.prop.lighting = False + self._actors["occlusion_surf"] = occl_act + elif surf_map["map_kind"] == "eeg": + head_mesh = surf_map["mesh"] + head_mesh._actor.prop.culling = "back" + self._actors["head"] = head_mesh._actor + + show_meg = (self._ch_type is None or self._ch_type == "meg") and any( + [m["kind"] == "meg" for m in self._surf_maps] + ) + show_eeg = (self._ch_type is None or self._ch_type == "eeg") and any( + [m["kind"] == "eeg" for m in self._surf_maps] + ) + meg_picks = pick_types(self._evoked.info, meg=show_meg, ref_meg=False) + eeg_picks = pick_types(self._evoked.info, meg=False, eeg=show_eeg) + picks = np.concatenate((meg_picks, eeg_picks)) + self._ch_names = [self._evoked.ch_names[i] for i in picks] + + for m in self._surf_maps: + if m["kind"] == "eeg": + head_surf = m["surf"] + break + else: + self._actors["head"], _, head_surf = _plot_head_surface( + renderer=fig_ef._renderer, + head="head", + subject=self._subject, + subjects_dir=self._subjects_dir, + bem=self._bem, + coord_frame="mri", + to_cf_t=self._to_cf_t, + alpha=0.2, + ) + self._actors["head"].prop.culling = "back" + + if show_sensors: + sensors = _plot_sensors_3d( + renderer=fig_ef._renderer, + info=self._evoked.info, + to_cf_t=self._to_cf_t, + picks=picks, + meg=["sensors"] if show_meg else False, + eeg=["original"] if show_eeg else False, + fnirs=False, + warn_meg=False, + head_surf=head_surf, + units="m", + sensor_alpha=dict(meg=0.1, eeg=1.0), + orient_glyphs=False, + scale_by_distance=False, + project_points=False, + surf=None, + check_inside=None, + nearest=None, + sensor_colors=dict( + meg=["gray" for _ in meg_picks], + eeg=["white" for _ in eeg_picks], + ), + ) + self._actors["sensors"] = sum(sensors.values(), []) + + # Adjust camera + fig_ef._renderer.set_camera( + azimuth=180, elevation=90, roll=90, distance=0.55, focalpoint=[0, 0, 0.03] + ) + + subscribe(fig_ef, "time_change", self._on_time_change) + subscribe(fig_ef, "channels_select", self._on_channels_select) + self._fig = fig_ef + + def _configure_dock(self): + """Configure the left and right dock areas of the GUI.""" + r = self._renderer + + # Toggle buttons for various meshes + layout = r._dock_add_group_box("Meshes", collapse=True) + + @_auto_weakref + def _toggle_mesh(_, name, show=None): + self.toggle_mesh(name, show=show) + + for actor_name in self._actors: + if actor_name == "occlusion_surf": + continue + r._dock_add_check_box( + name=actor_name, + value=True, + callback=partial(_toggle_mesh, name=actor_name), + layout=layout, + ) + + # Right dock + r._dock_initialize(name="Dipole fitting", area="right") + r._dock_add_button("Sensor data", self._on_sensor_data) + r._dock_add_button("Fit dipole", self.fit_dipole) + methods = ["Multi dipole (MNE)", "Single dipole"] + + @_auto_weakref + def _on_select_method(method): + self._on_select_method(method) + + r._dock_add_combo_box( + "Dipole model", + value="Multi dipole (MNE)", + rng=methods, + callback=_on_select_method, + ) + self._dipole_box = r._dock_add_group_box(name="Dipoles", collapse=False) + + @_auto_weakref + def _save(fname): + return self.save(fname) + + self._save_button = r._dock_add_file_button( + name="save_dipoles", + desc="Save dipoles", + save=True, + func=_save, + tooltip="Save the dipoles to disk", + filter_="Dipole files (*.dip *.bdip)", + initial_directory=".", + ) + self._save_button.set_enabled(False) + r._dock_add_stretch() + + def toggle_mesh(self, name, show=None): + """Toggle a mesh on or off. + + Parameters + ---------- + name : str + Name of the mesh to toggle. + show : bool | None + Whether to show the mesh. If None, the visibility of the mesh is toggled. + """ + _check_option("name", name, self._actors.keys()) + actors = self._actors[name] + # self._actors[name] is sometimes a list and sometimes not. Make it + # always be a list to simplify the code. + if not isinstance(actors, list): + actors = [actors] + if show is None: + show = not actors[0].GetVisibility() + for act in actors: + act.SetVisibility(show) + self._renderer._update() + + def set_time(self, time): + """Set the time point currently shown in the GUI. + + This is the programmatic equivalent of dragging the time slider, and is also the + time at which :meth:`fit_dipole` will fit a dipole. + + Parameters + ---------- + time : float + The time to show, in seconds. Values outside the time range of the evoked + data are clipped to the nearest valid time. + """ + publish(self._fig, TimeChange(time=float(time))) + + def _on_time_change(self, event): + new_time = np.clip(event.time, self._evoked.times[0], self._evoked.times[-1]) + self._current_time = new_time + print("gui time change to", new_time) + if self._time_line is not None: + self._time_line.set_xdata([new_time]) + self._renderer._mplcanvas.update_plot() + self._update_arrows() + + # TODO: Need to expose a public method for opening the sensor-data window and for + # programmatically selecting the channels to fit dipoles to. + def _on_sensor_data(self): + """Show sensor data and allow sensor selection.""" + if self._fig_sensors is not None: + return + fig = self._evoked.plot_topo(select=True) + fig.canvas.mpl_connect("close_event", self._on_sensor_data_close) + link(self._fig, fig, recursive=True) + self._fig_sensors = fig + + def _on_sensor_data_close(self, event): + """Handle closing of the sensor selection window.""" + publish(self._fig, ChannelsSelect(ch_names=[])) + self._fig_sensors = None + + def _on_channels_select(self, event): + """Color selected sensor meshes.""" + selected_channels = set(event.ch_names) + if "sensors" in self._actors: + # Possibly multiple sensor types. + for actor in self._actors["sensors"]: + cloud = actor.GetMapper().GetInput() + selected_idx = np.isin( + cloud.field_data["ch_names"], list(selected_channels) + ) + colors = cloud.point_data["colors"] + colors[selected_idx] = [0, 255, 0, 100] + colors[~selected_idx] = [0, 0, 0, 10] + cloud.point_data["colors"] = colors + self._renderer._update() + + def fit_dipole(self): + """Fit a single dipole and add it to the model. + + This is the programmatic equivalent of pressing the "Fit dipole" button. The + dipole is fitted at the time currently shown in the GUI (see :meth:`set_time`), + using the sensors that are currently selected in the sensor data window (or all + sensors when no selection is active). The newly fitted dipole is appended to the + :attr:`dipoles` attribute. + """ + evoked_picked = self._evoked.copy() + cov_picked = self._cov.copy() + if self._fig_sensors is not None: + picks = self._fig_sensors.lasso.selection + if len(picks) > 0: + evoked_picked = evoked_picked.pick(picks) + evoked_picked.info.normalize_proj() + cov_picked = cov_picked.pick_channels(picks, ordered=False) + cov_picked["projs"] = evoked_picked.info["projs"] + evoked_picked.crop(self._current_time, self._current_time) + + dip = fit_dipole( + evoked_picked, + cov_picked, + self._bem, + trans=self._head_mri_t, + rank=self._rank, + n_jobs=self._n_jobs, + verbose=False, + )[0] + + self.add_dipole(dip) + + def add_dipole(self, dipole, name=None): + """Add a dipole (or multiple dipoles) to the GUI. + + Parameters + ---------- + dipole : Dipole + The dipole to add. If the ``Dipole`` object defines multiple dipoles, they + will all be added. + name : str | list of str | None + The name of the dipole. When the ``Dipole`` object defines multiple dipoles, + this should be a list containing the name for each dipole. When ``None``, + the ``.name`` attribute of the ``Dipole`` object itself will be used. + """ + _validate_type(name, (str, list, None), "name") + if isinstance(name, str): + names = [name] + elif name is None: + # Try to obtain names from `dipole.name`. When multiple dipoles are saved, + # the names are concatenated with `;` marks. + if dipole.name is None: + names = [None] * len(dipole) + elif len(dipole.name.split(";")) == len(dipole): + names = dipole.name.split(";") + else: + names = [dipole.name] * len(dipole) + else: + names = name + if len(names) != len(dipole): + raise ValueError( + f"Number of names ({len(names)}) does not match the number of dipoles " + f"({len(dipole)})." + ) + + # Ensure orientations are unit vectors. Due to rounding issues this is sometimes + # not the case. + dipole._ori /= np.linalg.norm(dipole._ori, axis=1, keepdims=True) + + @_auto_weakref + def _on_dipole_toggle(active, dip_num): + return self._on_dipole_toggle(active, dip_num) + + @_auto_weakref + def _on_dipole_set_name(name, dip_num): + return self._on_dipole_set_name(name, dip_num) + + @_auto_weakref + def _on_dipole_toggle_fix_orientation(fix, dip_num): + return self._on_dipole_toggle_fix_orientation(fix, dip_num) + + @_auto_weakref + def _on_dipole_delete(dip_num): + return self._on_dipole_delete(dip_num) + + new_dipoles = list() + for dip, name in zip(dipole, names): + # Coordinates needed to draw the big arrow on the helmet. + helmet_coords, helmet_pos = self._get_helmet_coords(dip) + + # Collect all relevant information on the dipole in a dict. + colors = _get_color_list() + if len(self._dipoles) == 0: + dip_num = 0 + else: + dip_num = max(self._dipoles.keys()) + 1 + if name is None: + dip.name = f"dip{dip_num}" + else: + dip.name = name + dip_color = colors[dip_num % len(colors)] + if helmet_coords is not None: + arrow_mesh = pyvista.PolyData(*_arrow_mesh()) + else: + arrow_mesh = None + dipole_dict = dict( + active=True, + brain_arrow_actor=None, + helmet_arrow_actor=None, + arrow_mesh=arrow_mesh, + color=dip_color, + dip=dip, + fix_ori=True, + fix_position=True, + helmet_coords=helmet_coords, + helmet_pos=helmet_pos, + num=dip_num, + # fit_time=self._current_time, + ) + self._dipoles[dip_num] = dipole_dict + + # Add a row to the dipole list + r = self._renderer + hlayout = r._dock_add_layout(vertical=False) + widgets = [] + widgets.append( + r._dock_add_check_box( + name="", + value=True, + callback=partial(_on_dipole_toggle, dip_num=dip_num), + layout=hlayout, + ) + ) + widgets.append( + r._dock_add_text( + name=dip.name, + value=dip.name, + placeholder="name", + callback=partial(_on_dipole_set_name, dip_num=dip_num), + layout=hlayout, + ) + ) + widgets.append( + r._dock_add_check_box( + name="Fix ori", + value=True, + callback=partial( + _on_dipole_toggle_fix_orientation, dip_num=dip_num + ), + layout=hlayout, + ) + ) + widgets.append( + r._dock_add_button( + name="", + icon="clear", + callback=partial(_on_dipole_delete, dip_num=dip_num), + layout=hlayout, + ) + ) + dipole_dict["widgets"] = widgets + r._layout_add_widget(self._dipole_box, hlayout) + new_dipoles.append(dipole_dict) + + # Show the dipoles and arrows in the 3D view. Only do this after + # `_fit_timecourses` so that they have the correct size straight away. + self._fit_timecourses() + for dipole_dict in new_dipoles: + dip = dipole_dict["dip"] + dipole_dict["brain_arrow_actor"] = self._renderer.plotter.add_arrows( + apply_trans(self._head_mri_t, dip.pos[0]), + apply_trans(self._head_mri_t, dip.ori[0]), + color=dipole_dict["color"], + mag=0.05, + ) + if dipole_dict["arrow_mesh"] is not None: + dipole_dict["helmet_arrow_actor"] = self._renderer.plotter.add_mesh( + dipole_dict["arrow_mesh"], + color=dipole_dict["color"], + culling="front", + ) + self._update_arrows() + + def _get_helmet_coords(self, dip): + """Compute the coordinate system used for drawing the big arrows on the helmet. + + In this coordinate system, Z is normal to the helmet surface, and XY + are tangential to the helmet surface. + """ + if "helmet" not in self._actors: + return None, None + + # Get the closest vertex (=point) of the helmet mesh + dip_pos = apply_trans(self._head_mri_t, dip.pos[0]) + points = self._helmet_surf["rr"] + normals = self._helmet_surf["nn"] + distances = ((points - dip_pos) * normals).sum(axis=1) + closest_point = np.argmin(distances) + + # Compute the position of the projected dipole on the helmet + norm = normals[closest_point] + helmet_pos = dip_pos + (distances[closest_point] + 0.003) * norm + + # Create a coordinate system where X and Y are tangential to the helmet + helmet_coords = _normal_orth(norm) + + return helmet_coords, helmet_pos + + def _fit_timecourses(self): + """Compute (or re-compute) dipole timecourses. + + Called whenever something changes to the multi-dipole situation, i.e. a dipole + is added, removed, (de-)activated or the "Fix pos" box is toggled. + """ + self._save_button.set_enabled(len(self.dipoles) > 0) + active_dips = [d for d in self._dipoles.values() if d["active"]] + if len(active_dips) == 0: + return + + if self._multi_dipole_method == "Multi dipole (MNE)": + # TODO: When two active dipoles have (nearly) identical positions, they + # collapse to a single point in the discrete source space below, which + # errors out. Ideal behavior unclear: merge them, or error informatively? + this_src = setup_volume_source_space( + "sample", + pos=dict( + rr=apply_trans( + self._head_mri_t, + np.vstack([d["dip"].pos[0] for d in active_dips]), + ), + nn=apply_trans( + self._head_mri_t, + np.vstack([d["dip"].ori[0] for d in active_dips]), + ), + ), + ) + this_fwd = self.fwd.compute(this_src) + this_fwd = convert_forward_solution(this_fwd, surf_ori=False) + + inv = make_inverse_operator( + self._evoked.info, + # fwd, + this_fwd, + self._cov, + fixed=False, + loose=1.0, + depth=0, + rank=self._rank, + ) + stc = apply_inverse( + self._evoked, + inv, + method="MNE", + lambda2=1e-6, + pick_ori="vector", + ) + + timecourses = stc.magnitude().data + orientations = (stc.data / timecourses[:, np.newaxis, :]).transpose(0, 2, 1) + fixed_timecourses = stc.project( + np.array([dip["dip"].ori[0] for dip in active_dips]) + )[0].data + + for i, dip in enumerate(active_dips): + if dip["fix_ori"]: + dip["timecourse"] = fixed_timecourses[i] + dip["orientation"] = dip["dip"].ori.repeat(len(stc.times), axis=0) + else: + dip["timecourse"] = timecourses[i] + dip["orientation"] = orientations[i] + else: + assert self._multi_dipole_method == "Single dipole" # only other option + for dip in active_dips: + dip_with_timecourse, _ = fit_dipole( + self._evoked, + self._cov, + self._bem, + pos=dip["dip"].pos[0], # position is always fixed + ori=dip["dip"].ori[0] if dip["fix_ori"] else None, + trans=self._head_mri_t, + rank=self._rank, + n_jobs=self._n_jobs, + verbose=True, + ) + if dip["fix_ori"]: + dip["timecourse"] = dip_with_timecourse.data[0] + dip["orientation"] = dip["dip"].ori.repeat( + len(dip_with_timecourse.times), axis=0 + ) + else: + dip["timecourse"] = dip_with_timecourse.amplitude + dip["orientation"] = dip_with_timecourse.ori + + # Update matplotlib canvas at the bottom of the window + canvas = self._setup_mplcanvas() + ymin, ymax = 0, 0 + for dip in active_dips: + if "line_artist" in dip: + dip["line_artist"].set_ydata(dip["timecourse"]) + else: + dip["line_artist"] = canvas.plot( + self._evoked.times, + dip["timecourse"], + label=dip["dip"].name, + color=dip["color"], + ) + ymin = min(ymin, 1.1 * dip["timecourse"].min()) + ymax = max(ymax, 1.1 * dip["timecourse"].max()) + canvas.axes.set_ylim(ymin, ymax) + canvas.update_plot() + self._update_arrows() + + @verbose + def save(self, fname, verbose=None): + """Save the fitted dipoles to a file. + + Parameters + ---------- + fname : path-like + The name of the file. Should end in ``'.dip'`` to save in plain text format, + or in ``'.bdip'`` to save in binary format. + %(verbose)s + """ + if len(self.dipoles) == 0: + logger.info("No dipoles to save.") + return + + logger.info(f"Saving dipoles as: {fname}") + fname = Path(fname) + + # Pack the dipoles into a single mne.Dipole object. + if all(d.khi2 is not None for d in self.dipoles): + khi2 = np.array([d.khi2[0] for d in self.dipoles]) + else: + khi2 = None + + if all(d.nfree is not None for d in self.dipoles): + nfree = np.array([d.nfree[0] for d in self.dipoles]) + else: + nfree = None + + dip = Dipole( + times=np.array([d.times[0] for d in self.dipoles]), + pos=np.array([d.pos[0] for d in self.dipoles]), + amplitude=np.array([d.amplitude[0] for d in self.dipoles]), + ori=np.array([d.ori[0] for d in self.dipoles]), + gof=np.array([d.gof[0] for d in self.dipoles]), + khi2=khi2, + nfree=nfree, + conf={ + key: np.array([d.conf[key][0] for d in self.dipoles]) + for key in self.dipoles[0].conf.keys() + }, + name=";".join(d.name if hasattr(d, "name") else "" for d in self.dipoles), + ) + dip.save(fname, overwrite=True, verbose=verbose) + + def _update_arrows(self): + """Update the arrows to have the correct size and orientation.""" + active_dips = [d for d in self._dipoles.values() if d["active"]] + if len(active_dips) == 0: + return + orientations = [dip["orientation"] for dip in active_dips] + timecourses = [dip["timecourse"] for dip in active_dips] + arrow_scaling = 0.05 / np.max(np.abs(timecourses)) + for dip, ori, timecourse in zip(active_dips, orientations, timecourses): + helmet_coords = dip["helmet_coords"] + if helmet_coords is None: + continue + + dip_ori = apply_trans( + self._head_mri_t, + [np.interp(self._current_time, self._evoked.times, o) for o in ori.T], + ) + dip_moment = np.interp(self._current_time, self._evoked.times, timecourse) + arrow_size = dip_moment * arrow_scaling + arrow_mesh = dip["arrow_mesh"] + + # Project the orientation of the dipole tangential to the helmet + dip_ori_tan = helmet_coords[:2] @ dip_ori @ helmet_coords[:2] + + # Rotate the coordinate system such that Y lies along the dipole + # orientation, now we have our desired coordinate system for the + # arrows. + arrow_coords = np.array( + [np.cross(dip_ori_tan, helmet_coords[2]), dip_ori_tan, helmet_coords[2]] + ) + arrow_coords /= np.linalg.norm(arrow_coords, axis=1, keepdims=True) + + # Update the arrow mesh to point in the right directions + arrow_mesh.points = (_arrow_mesh()[0] * arrow_size) @ arrow_coords + arrow_mesh.points += dip["helmet_pos"] + self._renderer._update() + + # TODO: Need to expose a public method for setting the multi-dipole method + def _on_select_method(self, method): + """Select the method to use for multi-dipole timecourse fitting.""" + self._multi_dipole_method = method + self._fit_timecourses() + + # TODO: Need to expose public methods for toggling, renaming, (un)fixing the + # orientation of, and deleting a dipole (probably addressed by name or index). + def _on_dipole_toggle(self, active, dip_num): + """Toggle a dipole on or off.""" + dipole = self._dipoles[dip_num] + active = bool(active) + dipole["active"] = active + dipole["line_artist"].set_visible(active) + # Labels starting with "_" are hidden from the legend. + dipole["line_artist"].set_label(("" if active else "_") + dipole["dip"].name) + dipole["brain_arrow_actor"].visibility = active + dipole["helmet_arrow_actor"].visibility = active + self._fit_timecourses() + self._renderer._update() + self._renderer._mplcanvas.update_plot() + + def _on_dipole_set_name(self, name, dip_num): + """Set the name of a dipole.""" + self._dipoles[dip_num]["dip"].name = name + self._dipoles[dip_num]["line_artist"].set_label(name) + self._renderer._mplcanvas.update_plot() + + def _on_dipole_toggle_fix_orientation(self, fix, dip_num): + """Fix dipole orientation when fitting timecourse.""" + self._dipoles[dip_num]["fix_ori"] = bool(fix) + self._fit_timecourses() + + def _on_dipole_delete(self, dip_num): + """Delete previously fitted dipole.""" + dipole = self._dipoles[dip_num] + dipole["line_artist"].remove() + dipole["brain_arrow_actor"].visibility = False + if dipole["helmet_arrow_actor"] is not None: # no helmet arrow for EEG + dipole["helmet_arrow_actor"].visibility = False + for widget in dipole["widgets"]: + widget.hide() + del self._dipoles[dip_num] + self._fit_timecourses() + self._renderer._update() + self._renderer._mplcanvas.update_plot() + + def _setup_mplcanvas(self): + """Configure the matplotlib canvas at the bottom of the window.""" + if self._renderer._mplcanvas is None: + self._renderer._mplcanvas = self._renderer._window_get_mplcanvas( + self._fig, 0.3, False, False + ) + self._renderer._window_adjust_mplcanvas_layout() + if self._time_line is None: + self._time_line = self._renderer._mplcanvas.plot_time_line( + self._current_time, + label="time", + color="black", + ) + return self._renderer._mplcanvas + + def close(self): + """Close the dipole fitting GUI.""" + if self._renderer is not None: + try: + self._renderer.close() + except AttributeError: # maybe already closed + pass + + +def _arrow_mesh(): + """Obtain a mesh of an arrow.""" + vertices = np.array( + [ + [0.0, 1.0, 0.0], + [0.3, 0.7, 0.0], + [0.1, 0.7, 0.0], + [0.1, -1.0, 0.0], + [-0.1, -1.0, 0.0], + [-0.1, 0.7, 0.0], + [-0.3, 0.7, 0.0], + ] + ) + faces = np.array([[7, 0, 1, 2, 3, 4, 5, 6]]) + return vertices, faces diff --git a/mne/gui/_gui.py b/mne/gui/_gui.py index b8898d8b7c2..6f3b5474b83 100644 --- a/mne/gui/_gui.py +++ b/mne/gui/_gui.py @@ -165,16 +165,240 @@ def coregistration( ) +@verbose +def dipolefit( + evoked, + *, + baseline=None, + cov=None, + bem=None, + initial_time=None, + trans=None, + stc=None, + subject=None, + subjects_dir=None, + surf_maps=None, + rank="info", + show_density=True, + ch_type=None, + show_sensors=True, + n_jobs=None, + show=True, + block=False, + verbose=None, +): + """GUI for interactive dipole fitting, inspired by MEGIN's XFit program. + + Parameters + ---------- + evoked : instance of Evoked | path-like | None + Evoked data to show fieldmap of and fit dipoles to. + %(baseline_evoked)s + cov : instance of Covariance | path-like | "baseline" | None + Noise covariance matrix. If ``None``, an ad-hoc covariance matrix is used with + default values for the diagonal elements (see Notes). If ``"baseline"``, the + diagonal elements is estimated from the baseline period of the evoked data. + bem : instance of ConductorModel | path-like | None + Boundary element model to use in forward calculations, or a path to the BEM + solution file (``"-bem-sol.fif"``) to read it from. If ``None``, a spherical + model is used. + initial_time : float | None + Initial time point to show. If ``None``, the time point of the maximum field + strength is used. + trans : instance of Transform | path-like | None + The transformation from head coordinates to MRI coordinates. If ``None``, + the identity matrix is used and everything will be done in head coordinates. + stc : instance of SourceEstimate | path-like | None + An optional distributed source estimate to show alongside the fieldmap. The time + samples need to match those of the evoked data. + subject : str | None + The subject name. If ``None``, no MRI data is shown. + %(subjects_dir)s + surf_maps : list | None + The surface mapping information obtained with make_field_map. If ``None``, one + will be generated based on the given data. + %(rank)s + show_density : bool + Whether to show the density of the fieldmap. + ch_type : "meg" | "eeg" | None + Type of channels to use for the dipole fitting. By default (``None``) both MEG + and EEG channels will be used. + show_sensors : bool + Whether to show the sensors in the 3D view. + %(n_jobs)s + show : bool + Show the GUI if True. + block : bool + Whether to halt program execution until the figure is closed. + %(verbose)s + + Returns + ------- + fitter : instance of DipoleFitUI + The dipole fitting GUI. The ``.dipoles`` attribute contains the fitted dipoles. + + Notes + ----- + .. versionadded:: 1.12 + + When using ``cov=None`` the default noise values are 5 fT/cm, 20 fT, and 0.2 µV for + gradiometers, magnetometers, and EEG channels respectively. + + Here is an incomplete comparison between the features of the MEGIN XFit™ 5.5.18 + software and the MNE interactive dipole fitting GUI: + + .. table:: + :widths: auto + + +-----------------------------------------------------------------------------+-----+------+ + | Feature | MNE | Xfit | + +=============================================================================+=====+======+ + | | + +-----------------------------------------------------------------------------+-----+------+ + | **Head model** | + +-----------------------------------------------------------------------------+-----+------+ + | Use spherical head model | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Use BEM head model based on MRI | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Use ad-hoc covariance matrix | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Estimate covariance from baseline period | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | **Dipole fitting** | + +-----------------------------------------------------------------------------+-----+------+ + | Fit a dipole at the current time | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Fit a dipole by averaging the signal over a time range | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Fit dipoles on a subset of sensors, selected from the sensor view | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Fit multiple dipoles and construct a multi-dipole model | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Toggle individual dipoles on or off in the multi-dipole model | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Give names to dipoles | ✓ | | + +-----------------------------------------------------------------------------+-----+------+ + | Save dipoles to .dip or .bdip file | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | View source timecourses | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | View total variance explained by the dipole model | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | View detailed dipole information (coordinates, goodness of fit, etc.) | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Toggle dipoles to have a free or fixed orientation | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Project-out signals from dipoles currently in the model | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | **3D view** | + +-----------------------------------------------------------------------------+-----+------+ + | View magnetic field patterns for Evokeds | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Show head surface in relation to the MEG helmet | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Show brain surfaces inside the head surface | ✓ | | + +-----------------------------------------------------------------------------+-----+------+ + | Show MNE source estimate to guide dipole fitting | ✓ | | + +-----------------------------------------------------------------------------+-----+------+ + | Show location of the dipoles in the 3D view | ✓ | | + +-----------------------------------------------------------------------------+-----+------+ + | Show dipole projected to the helmet (big arrows) | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Show currently selected sensors in the 3D view | ✓ | | + +-----------------------------------------------------------------------------+-----+------+ + | Only show magnetic field patterns as seen by the currently selected sensors | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | **Sensor view** | + +-----------------------------------------------------------------------------+-----+------+ + | View sensor-level timecourses | ✓ | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Switch layouts (all, grads, mags, eeg) in the sensor display | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + | Overlay reconstructed sensor time courses from the dipole model | | ✓ | + +-----------------------------------------------------------------------------+-----+------+ + """ # noqa E501 + from ..viz.backends.renderer import MNE_3D_BACKEND_TESTING + from ._dipolefit import DipoleFitUI + + if MNE_3D_BACKEND_TESTING: + show = block = False + + return DipoleFitUI( + evoked=evoked, + baseline=baseline, + cov=cov, + bem=bem, + initial_time=initial_time, + trans=trans, + stc=stc, + subject=subject, + subjects_dir=subjects_dir, + surf_maps=surf_maps, + rank=rank, + show_density=show_density, + ch_type=ch_type, + show_sensors=show_sensors, + n_jobs=n_jobs, + show=show, + block=block, + verbose=verbose, + ) + + +def _gui_closed(gui): + """Check whether a GUI's 3D renderer has already been closed.""" + if not hasattr(gui, "_renderer"): # nothing to close + return False + try: + plotter = gui._renderer.plotter + except Exception: + return True + return bool(getattr(plotter, "_closed", False)) + + +def _close_gui(gui): + """Close a GUI, tolerating GUIs that are already (partially) closed.""" + try: # for compatibility with both GUIs, will be refactored + gui._renderer.close() # TODO should be triggered by close + except Exception: + pass + gui.close() + + class _GUIScraper: - """Scrape GUI outputs.""" + """Scrape GUI outputs. + + By default a GUI is scraped once and then closed, so each GUI shows up as a single + image in the rendered example. An example can instead add a file-level + + .. code-block:: python + + # sphinx_gallery_preserve_gui = True + + comment, in which case the GUI is scraped after *every* code block and left open, so + that successive code blocks can keep operating on it. The scraper then holds the + only strong reference to those GUIs, and ``close_preserved`` (called from the doc + build's ``reset_modules``) closes them once the example is done. + """ + + def __init__(self): + self._preserved_guis = list() def __repr__(self): return "" + def close_preserved(self): + """Close (and forget about) all GUIs preserved across code blocks.""" + guis, self._preserved_guis = self._preserved_guis, list() + for gui in guis: + _close_gui(gui) + def __call__(self, block, block_vars, gallery_conf): from ._coreg import CoregistrationUI + from ._dipolefit import DipoleFitUI - gui_classes = (CoregistrationUI,) + gui_classes = (CoregistrationUI, DipoleFitUI) try: from mne_gui_addons._ieeg_locate import IntracranialElectrodeLocator except Exception: @@ -184,11 +408,18 @@ def __call__(self, block, block_vars, gallery_conf): from qtpy import QtGui from sphinx_gallery.scrapers import figure_rst + preserve = bool((block_vars.get("file_conf") or {}).get("preserve_gui", False)) for gui in block_vars["example_globals"].values(): if ( isinstance(gui, gui_classes) - and not getattr(gui, "_scraped", False) and gallery_conf["builder_name"] == "html" + and ( + # A preserved GUI is scraped for every code block; any other one + # only the first time we see it. + not _gui_closed(gui) + if preserve + else not getattr(gui, "_scraped", False) + ) ): gui._scraped = True # monkey-patch but it's easy enough img_fname = next(block_vars["image_path_iterator"]) @@ -213,10 +444,19 @@ def __call__(self, block, block_vars, gallery_conf): ) # https://doc.qt.io/qt-5/qpixmap.html#save pixmap.save(img_fname) - try: # for compatibility with both GUIs, will be refactored - gui._renderer.close() # TODO should be triggered by close - except Exception: - pass - gui.close() + if preserve: + # Keep it open (and alive) so the next code block can use it. + if not any(gui is known for known in self._preserved_guis): + self._preserved_guis.append(gui) + if hasattr(gui, "_renderer"): + # The PyVista scraper runs after us and screenshots *and then + # closes* every plotter it knows about, which would take our GUI + # down with it. Deregister ours from it, just like closing a + # figure does (see _pyvista._close_3d_figure). + from ..viz.backends._pyvista import _ALL_PLOTTERS + + _ALL_PLOTTERS.pop(plotter._id_name, None) + else: + _close_gui(gui) return figure_rst([img_fname], gallery_conf["src_dir"], "GUI") return "" diff --git a/mne/gui/tests/test_dipolefit.py b/mne/gui/tests/test_dipolefit.py new file mode 100644 index 00000000000..a3bfc8a8e9b --- /dev/null +++ b/mne/gui/tests/test_dipolefit.py @@ -0,0 +1,523 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_equal + +import mne +from mne.channels import read_vectorview_selection +from mne.datasets import testing +from mne.viz import ui_events +from mne.viz.utils import _get_color_list + +data_path = testing.data_path(download=False) +subjects_dir = data_path / "subjects" +fname_dip = data_path / "MEG" / "sample" / "sample_audvis_trunc_set1.dip" +fname_evokeds = data_path / "MEG" / "sample" / "sample_audvis_trunc-ave.fif" +fname_trans = data_path / "MEG" / "sample" / "sample_audvis_trunc-trans.fif" +fname_cov = data_path / "MEG" / "sample" / "sample_audvis_trunc-cov.fif" +fname_stc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg" +fname_bem_sol = subjects_dir / "sample" / "bem" / "sample-320-320-320-bem-sol.fif" + + +@pytest.fixture(scope="module") +def _sample_evoked(): + """Read the evoked data (module scoped, do not modify: use ``sample_evoked``).""" + return mne.read_evokeds(fname_evokeds, condition=0) + + +@pytest.fixture +def sample_evoked(_sample_evoked): + """Get the evoked data used throughout these tests.""" + return _sample_evoked.copy() + + +@pytest.fixture(scope="module") +def surf_maps_meg(_sample_evoked): + """Compute the MEG-only field map (as ``dipolefit`` would without a trans).""" + return mne.make_field_map(_sample_evoked, trans=None, origin="auto", verbose=False) + + +@pytest.fixture(scope="module") +def surf_maps_eeg_meg(_sample_evoked): + """Compute both the EEG and MEG field maps (needs a head<->MRI transform).""" + return mne.make_field_map( + _sample_evoked, + trans=fname_trans, + origin="auto", + subject="sample", + subjects_dir=subjects_dir, + verbose=False, + ) + + +def _selected_sensors(g): + names = [] + for actor in g._actors["sensors"]: + cloud = actor.GetMapper().GetInput() + # color is hardcoded for now, so changeable + green = (cloud.point_data["colors"] == [0, 255, 0, 100]).all(axis=1) + names.extend(cloud.field_data["ch_names"][green]) + return sorted(names) + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_gui_basic( + tmp_path, sample_evoked, surf_maps_meg, renderer_interactive_pyvistaqt +): + """Test basic functionality of the dipole fitting GUI.""" + from mne.gui import dipolefit + + # Test basic interface elements. + evoked = sample_evoked + data_before = evoked.data.copy() + g = dipolefit(evoked, baseline=(None, 0), surf_maps=surf_maps_meg) + + assert evoked.comment == "Left Auditory" # MNE-Sample data should be loaded + assert_allclose(evoked.data, data_before, atol=0) # input is not modified + assert g._evoked.baseline == (evoked.times[0], 0) # baseline applied to a copy + assert g._current_time == evoked.times[84] # time of max GFP + + # The sensors consist of multiple actors, toggling them should affect all of them. + g.toggle_mesh("sensors", show=False) + assert not g._actors["sensors"][0].GetVisibility() + g.toggle_mesh("sensors") # show=None toggles the current visibility + assert g._actors["sensors"][0].GetVisibility() + + # Test fitting a single dipole. + assert len(g._dipoles) == len(g.dipoles) == 0 + g.fit_dipole() + assert len(g._dipoles) == len(g.dipoles) == 1 + dip = g.dipoles[0] + assert dip.name == "Left Auditory" + assert len(dip.times) == 1 + assert_equal(dip.times, g._current_time) + old_dip1_timecourse = g._dipoles[0]["timecourse"] + + # Check the position of the fitted dipole against the pre-computed dipole in the + # testing dataset. The pre-computed dipole only needs to give us the general area + # in which we expect the fitted dipole and does not have to be a perfect match. + ref_dip = mne.read_dipole(fname_dip) + ref_dip_t = ref_dip[[np.argmin(np.abs(ref_dip.times - g._current_time))]] + assert_allclose(dip.pos, ref_dip_t.pos, atol=0.012) # somewhat near the reference + + # Test fitting a second dipole with a subset of channels at a different time. + g._on_sensor_data() # open sensor selection window + g._on_sensor_data() # already open, so this is a no-op + picks = read_vectorview_selection("Left", info=evoked.info) + ui_events.publish(g._fig_sensors, ui_events.ChannelsSelect(picks)) + assert sorted(g._fig_sensors.lasso.selection) == sorted(picks) + assert _selected_sensors(g) == sorted(picks) + ui_events.publish(g._fig, ui_events.TimeChange(0.09)) # change time + assert g._current_time == 0.09 + g.fit_dipole() + assert len(g._dipoles) == len(g.dipoles) == 2 + dip2 = g.dipoles[1] + + # During tests, matplotlib does not open an actual window so we need to force the + # close event. + g._fig_sensors.canvas.callbacks.process("close_event", None) + assert _selected_sensors(g) == [] + + # The selected time of 0.09 is not actually in evoked.times, find the closest value + # that is (0.08990784...). That should be the time recorded in the dipole object. + closest_time = evoked.times[np.argmin(np.abs(evoked.times - g._current_time))] + assert dip2.times[0] == closest_time + + # Check that the general area of the second dipole is now in the left hemisphere. + ref_dip_t = ref_dip[[np.argmin(np.abs(ref_dip.times - g._current_time))]] + assert_allclose(dip2.pos, ref_dip_t.pos, atol=0.012) # somewhat near the reference + + # Adding the second dipole should have affected the timecourse of the first. + new_dip1_timecourse = g._dipoles[0]["timecourse"] + assert not np.allclose(old_dip1_timecourse, new_dip1_timecourse, atol=1e-10) + + # Test differences between the two dipoles + assert list(g._dipoles.keys()) == [0, 1] + dip1_dict, dip2_dict = g._dipoles.values() + assert dip1_dict["dip"] is dip + assert dip2_dict["dip"] is dip2 + assert dip1_dict["num"] == 0 + assert dip2_dict["num"] == 1 + assert dip1_dict["color"] == _get_color_list()[0] + assert dip2_dict["color"] == _get_color_list()[1] + + # Fitted dipoles have goodness-of-fit information that should be saved along. + fname = tmp_path / "fitted.dip" + g.save(fname) + assert mne.read_dipole(fname).khi2 is not None + + # Test changing dipole model + assert g._multi_dipole_method == "Multi dipole (MNE)" + old_timecourses = np.vstack((dip1_dict["timecourse"], dip2_dict["timecourse"])) + g._on_select_method("Single dipole") + new_timecourses = np.vstack((dip1_dict["timecourse"], dip2_dict["timecourse"])) + assert not np.allclose(old_timecourses, new_timecourses, atol=1e-10) + + g.close() + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_gui_dipole_controls( + sample_evoked, surf_maps_meg, renderer_interactive_pyvistaqt +): + """Test the controls for the dipoles in the dipole fitting GUI.""" + from mne.gui import dipolefit + + evoked = sample_evoked + g = dipolefit(evoked, surf_maps=surf_maps_meg, show_sensors=False) + + # Test toggling the visibility of the meshes. + assert list(g._actors.keys()) == ["helmet", "occlusion_surf", "head"] + g.toggle_mesh("helmet", show=True) + assert g._actors["helmet"].visibility + g.toggle_mesh("helmet") + assert not g._actors["helmet"].visibility + with pytest.raises(ValueError, match="Invalid value for the 'name' parameter"): + g.toggle_mesh("non existent") + + # Test toggling dipoles off and on. This is done through the GUI widgets, which are + # ordered: [active, name, fix orientation, delete]. + dip = mne.read_dipole(fname_dip)[[12, 15]] # 80ms and 90ms + g.add_dipole(dip, name=["rh", "lh"]) + dip1, dip2 = g._dipoles.values() + assert dip1["active"] and dip2["active"] + old_timecourses = np.vstack((dip1["timecourse"], dip2["timecourse"])) + dip2["widgets"][0].set_value(False) + assert not dip2["active"] + new_timecourses = np.vstack((dip1["timecourse"], dip2["timecourse"])) + assert not np.allclose(old_timecourses, new_timecourses, atol=1e-9) + + # With all dipoles disabled, there is nothing to fit and no arrows to update. + dip1["widgets"][0].set_value(False) + assert g.dipoles == [] + g.set_time(0.05) + assert g._current_time == 0.05 + + dip1["widgets"][0].set_value(True) + dip2["widgets"][0].set_value(True) + assert dip1["active"] and dip2["active"] + new_timecourses = np.vstack((dip1["timecourse"], dip2["timecourse"])) + assert np.allclose(old_timecourses, new_timecourses, atol=0) + + # Toggle fixed orientation off and on. + assert dip1["fix_ori"] and dip2["fix_ori"] + dip1["widgets"][2].set_value(False) + assert not dip1["fix_ori"] + new_timecourses = np.vstack((dip1["timecourse"], dip2["timecourse"])) + assert not np.allclose(old_timecourses, new_timecourses, atol=1e-9) + dip1["widgets"][2].set_value(True) + assert dip1["fix_ori"] + new_timecourses = np.vstack((dip1["timecourse"], dip2["timecourse"])) + assert np.allclose(old_timecourses, new_timecourses, atol=0) + + # Change the names of the dipoles. + dip1["widgets"][1].set_value("dipole1") + g._on_dipole_set_name("dipole2", dip2["num"]) + assert dip1["dip"].name == "dipole1" + assert dip2["dip"].name == "dipole2" + assert dip1["line_artist"].get_label() == "dipole1" # legend labels + assert dip2["line_artist"].get_label() == "dipole2" + + # Remove a dipole (through the "delete" button). + dip1["widgets"][3].set_value(None) + assert len(g.dipoles) == 1 + assert 1 in g._dipoles # dipole number should not change + assert list(g._dipoles.keys())[0] == 1 + assert list(g._dipoles.values())[0]["num"] == 1 + g.fit_dipole() + assert 2 in g._dipoles + assert list(g._dipoles.keys())[1] == 2 + assert list(g._dipoles.values())[1]["num"] == 2 # new dipole number + + # Fitting the timecourse of a single dipole, with a free orientation. + g._on_dipole_toggle(False, 2) # only leave a single dipole active + g._on_select_method("Single dipole") + assert dip2["fix_ori"] + assert_allclose(dip2["orientation"], dip2["dip"].ori.repeat(len(evoked.times), 0)) + g._on_dipole_toggle_fix_orientation(False, dip2["num"]) + assert not dip2["fix_ori"] + assert dip2["orientation"].shape == (len(evoked.times), 3) + assert not np.allclose(dip2["orientation"][0], dip2["orientation"][-1], atol=1e-9) + + g.close() + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_gui_save_load( + tmp_path, sample_evoked, renderer_interactive_pyvistaqt +): + """Test saving and loading dipoles in the dipole fitting GUI.""" + from mne.gui import dipolefit + + # Not passing `surf_maps` means they are computed on the fly. + g = dipolefit(sample_evoked, show_sensors=False) + dip = mne.read_dipole(fname_dip)[[12, 15]] # 80ms and 90ms + g.add_dipole(dip, name=["rh", "lh"]) + + g.save(tmp_path / "test.dip") + g.save(tmp_path / "test.bdip") + dip_from_file = mne.read_dipole(tmp_path / "test.dip") + g.add_dipole(dip_from_file) # names are taken from the ";" separated dip.name + g.add_dipole(mne.read_dipole(tmp_path / "test.bdip")) # bdip stores no names + assert len(g.dipoles) == 6 + assert [d.name for d in g.dipoles] == ["rh", "lh", "rh", "lh", "dip4", "dip5"] + for start in [0, 2, 4]: + assert_allclose( + np.vstack([d.pos for d in g.dipoles[start : start + 2]]), + dip_from_file.pos, + atol=0, + ) + + # A single dipole can be given a name directly. When the name of the `Dipole` object + # cannot be split into one name per dipole, it is used for all of them. + assert dip_from_file.name == "rh;lh" # cannot be split into a single name + g.add_dipole(dip_from_file[[0]], name="single") + g.add_dipole(dip_from_file[[1]]) + assert [d.name for d in g.dipoles[6:]] == ["single", "rh;lh"] + + with pytest.raises(ValueError, match="Number of names"): + g.add_dipole(dip_from_file, name=["too", "many", "names"]) + + g.close() + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_params( + tmp_path, sample_evoked, surf_maps_eeg_meg, renderer_interactive_pyvistaqt +): + """Test setting various parameters in the dipole fitting GUI.""" + from mne.gui import dipolefit + + # Test different type of covariance estimators. + evoked = sample_evoked + + g = dipolefit( + evoked, surf_maps=surf_maps_eeg_meg, cov=None, show_sensors=False + ) # ad-hoc + assert g._cov["diag"] + assert_allclose( # default ad-hoc variation for grads, mags and eeg + g._cov["data"][[0, 1, 2, 306]], [2.5e-25, 2.5e-25, 4e-28, 4e-14], atol=0 + ) + + # cov="baseline" needs baseline-corrected data (this evoked has baseline=None). + with pytest.raises(ValueError, match='cov="baseline" requires'): + dipolefit(evoked, surf_maps=surf_maps_eeg_meg, cov="baseline") + g = dipolefit( + evoked, + baseline=(None, 0), + surf_maps=surf_maps_eeg_meg, + cov="baseline", + show_sensors=False, + ) + assert_allclose( # compute var on baseline period + g._cov["data"][[0, 1, 2, 306]], + [3.5e-24, 3.5e-24, 3.0e-27, 2.3e-12], + rtol=0.1, + atol=0, + ) + + # The following tests are rolled into one call to `dipolefit` in order to save time. + # - Specify a channel type + # - Specify custom covariance. + # - Specify BEM model. + # - Specify an initial time. + cov = mne.read_cov(fname_cov) + bem = mne.make_sphere_model(r0=(0.0, 0.0, 0.04), verbose=False) + initial_time = 0.0123 + eeg_maps = [m for m in surf_maps_eeg_meg if m["kind"] == "eeg"] + g = dipolefit( + evoked, + ch_type="eeg", + surf_maps=eeg_maps, + cov=cov, + bem=bem, + initial_time=initial_time, + show_sensors=False, + ) + assert set(g._evoked.get_channel_types()) == {"eeg"} + assert_allclose(g._cov["data"], cov["data"], atol=0) + assert_equal(g._bem["r0"], bem["r0"]) + assert g._current_time == initial_time + + # Without an MEG helmet, no arrows are drawn on it. + g.add_dipole(mne.read_dipole(fname_dip)[[12]]) + (dipole,) = g._dipoles.values() + assert dipole["helmet_coords"] is None + assert dipole["arrow_mesh"] is None + assert dipole["helmet_arrow_actor"] is None + g._on_dipole_delete(dipole["num"]) + assert len(g.dipoles) == 0 + + # Without any dipoles, there is nothing to save. + g.save(tmp_path / "empty.dip") + assert not (tmp_path / "empty.dip").exists() + + g.close() + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_stc( + sample_evoked, surf_maps_eeg_meg, renderer_interactive_pyvistaqt +): + """Test showing a SourceEstimate underneath the fieldlines.""" + from mne.gui import dipolefit + + evoked = sample_evoked + + # By default, the STC file has different timestamps from the evoked. + with pytest.raises(ValueError, match="The time samples of the source estimate"): + dipolefit( + evoked, stc=fname_stc, surf_maps=surf_maps_eeg_meg, show_sensors=False + ) + + # Make the evoked timestamps line up with those of the STC. + stc = mne.read_source_estimate(fname_stc) + with pytest.warns(): + evoked = evoked.crop(0, 0.245).decimate(3) + with evoked.info._unlock(): + evoked.info["sfreq"] = 100 + evoked._set_times(stc.times) + + # A source estimate needs a transform to be shown in the correct place. + with pytest.raises(ValueError, match="`trans` cannot be `None`"): + dipolefit( + evoked, + stc=stc, + surf_maps=surf_maps_eeg_meg, + show_sensors=False, + baseline=(0, 0), + ) + + # Now it should work. Passing `bem` as a path loads the BEM solution file. + g = dipolefit( + evoked, + stc=stc, + trans=fname_trans, + subject="sample", + subjects_dir=subjects_dir, + surf_maps=surf_maps_eeg_meg, + bem=fname_bem_sol, + show_sensors=False, + baseline=(0, 0), + ) + assert isinstance(g._stc, mne.SourceEstimate) + assert not g._bem["is_sphere"] + assert "solution" in g._bem + g.close() + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_gui_scraper( + tmp_path, sample_evoked, surf_maps_meg, renderer_interactive_pyvistaqt +): + """Test the scraper for the dipole fitting GUI.""" + pytest.importorskip("sphinx_gallery") + from mne.gui import dipolefit + from mne.viz.backends._pyvista import _ALL_PLOTTERS + + (tmp_path / "_images").mkdir() + gallery_conf = dict(builder_name="html", src_dir=tmp_path) + scraper = mne.gui._GUIScraper() + + # By default a GUI is scraped once and then closed. + g = dipolefit(sample_evoked, surf_maps=surf_maps_meg, show_sensors=False) + img = tmp_path / "_images" / "temp.png" + block_vars = dict(example_globals=dict(gui=g), image_path_iterator=iter([str(img)])) + assert not getattr(g, "_scraped", False) + assert scraper(None, block_vars, gallery_conf) + assert img.is_file() + assert g._scraped + assert g._renderer.plotter._closed + assert scraper._preserved_guis == [] + assert scraper(None, block_vars, gallery_conf) == "" # only scraped once + + # With ``# sphinx_gallery_preserve_gui = True`` it is scraped for every code block + # and kept open until close_preserved() is called (from the doc build's + # reset_modules). + g = dipolefit(sample_evoked.copy(), surf_maps=surf_maps_meg, show_sensors=False) + assert g._renderer.plotter._id_name in _ALL_PLOTTERS + imgs = [tmp_path / "_images" / f"preserved{ii}.png" for ii in range(2)] + block_vars = dict( + example_globals=dict(gui=g), + image_path_iterator=iter([str(img) for img in imgs]), + file_conf=dict(preserve_gui=True), + ) + for img in imgs: + assert scraper(None, block_vars, gallery_conf) + assert img.is_file() + assert not g._renderer.plotter._closed + # deregistered from the PyVista scraper, which would otherwise screenshot the + # plotter a second time and then close it + assert g._renderer.plotter._id_name not in _ALL_PLOTTERS + assert scraper._preserved_guis == [g] + scraper.close_preserved() + assert scraper._preserved_guis == [] + assert g._renderer.plotter._closed + + +@pytest.mark.slowtest +@testing.requires_testing_data +def test_dipolefit_rapid_time_changes( + sample_evoked, surf_maps_eeg_meg, renderer_interactive_pyvistaqt +): + """Test that rapid time changes leave all linked views at the same time.""" + from qtpy.QtCore import QEvent, QObject + from qtpy.QtWidgets import QApplication + + from mne.gui import dipolefit + + # Same stc-aligned configuration as test_dipolefit_stc. + evoked = sample_evoked + stc = mne.read_source_estimate(fname_stc) + evoked = evoked.crop(0, 0.245).decimate(3, verbose="error") + with evoked.info._unlock(): + evoked.info["sfreq"] = 100 + evoked._set_times(stc.times) + g = dipolefit( + evoked, + stc=stc, + trans=fname_trans, + subject="sample", + subjects_dir=subjects_dir, + surf_maps=surf_maps_eeg_meg, + show_sensors=False, + baseline=(0, 0), + ) + + class _QueuedTimeChange(QEvent): + def __init__(self, time): + super().__init__(QEvent.Type.User) + self.time = time + + class _Publisher(QObject): + """Publish a TimeChange per Qt event, like rapid moves of the time slider.""" + + def customEvent(self, event): + # Publish on the brain figure so that the brain's handler (which + # processes pending Qt events, and thereby the next queued + # publication) runs before the other subscribers do. + ui_events.publish(g._stc_brain, ui_events.TimeChange(time=event.time)) + + # Queue up the publications as pending Qt events, then deliver them all, as + # happens when the time slider is scrolled faster than the views can redraw. + times = evoked.times[[10, 11, 12, 13]] + app = QApplication.instance() + publisher = _Publisher() + for time in times: + app.postEvent(publisher, _QueuedTimeChange(time)) + g._renderer._process_events() + + # Every view must end up at the last published time. + assert g._stc_brain._current_time == times[-1] # brain data + time line + assert g._fig._current_time == times[-1] # field lines + assert g._current_time == times[-1] # dipole arrows + g.close() diff --git a/mne/utils/misc.py b/mne/utils/misc.py index 9522c501a06..a25fc6b3e11 100644 --- a/mne/utils/misc.py +++ b/mne/utils/misc.py @@ -361,7 +361,6 @@ def _assert_no_instances(cls, when=""): from refleak.testing import assert_no_instances __tracebackhide__ = True - assert_no_instances(cls, when=when) diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py index 273a9f7e689..2d097d453a5 100644 --- a/mne/viz/_3d.py +++ b/mne/viz/_3d.py @@ -1342,7 +1342,7 @@ def _plot_hpi_coils( backface_culling=True, check_inside=check_inside, nearest=nearest, - ) + )[0] def _get_nearest(nearest, check_inside, project_to_trans, proj_rr): @@ -1413,7 +1413,7 @@ def _plot_glyphs( defaults = DEFAULTS["coreg"] n = len(loc) if n == 0: - return None + return None, None colors = np.array(np.broadcast_to(to_rgba_array(colors), (n, 4)), float) colors[:, 3] *= opacity scales = np.broadcast_to(np.asarray(scales, float).reshape(-1), (n,)) @@ -1445,7 +1445,7 @@ def _plot_glyphs( rots = np.array([_find_vector_rotation(x_axis, this_nn) for this_nn in nn]) quats = rot_to_quat(rots) rr, tris = renderer._glyph_template(kind, **template_kw) - actor, _ = renderer.instanced_mesh( + actor, cloud = renderer.instanced_mesh( rr=rr, tris=tris, positions=positions, @@ -1454,7 +1454,7 @@ def _plot_glyphs( scales=scales, backface_culling=backface_culling, ) - return actor + return actor, cloud @verbose @@ -1498,7 +1498,7 @@ def _plot_head_shape_points( backface_culling=True, check_inside=check_inside, nearest=nearest, - ) + )[0] def _plot_forward(renderer, fwd, fwd_trans, fwd_scale=1, scale=1.5e-3, alpha=1): @@ -1564,6 +1564,7 @@ def _plot_sensors_3d( actors = defaultdict(lambda: list()) locs = defaultdict(lambda: list()) + ch_names_all = defaultdict(lambda: list()) unit_scalar = 1 if units == "m" else 1e3 for ch_name, ch_coord in ch_pos.items(): ch_type = channel_type(info, info.ch_names.index(ch_name)) @@ -1597,14 +1598,19 @@ def _plot_sensors_3d( if ch_type == "eeg": if "original" in eeg: locs[ch_type].append(ch_coord) + ch_names_all[ch_type].append(ch_name) if "projected" in eeg: locs["eegp"].append(ch_coord) + ch_names_all["eegp"].append(ch_name) else: locs[ch_type].append(ch_coord) + ch_names_all[ch_type].append(ch_name) if ch_name in sources and "sources" in fnirs: locs["source"].append(sources[ch_name]) + ch_names_all["source"].append(ch_name) if ch_name in detectors and "detectors" in fnirs: locs["detector"].append(detectors[ch_name]) + ch_names_all["detector"].append(ch_name) # Plot these now if ch_name in sources and ch_name in detectors and "pairs" in fnirs: actor, _ = renderer.tube( # array of origin and dest points @@ -1665,7 +1671,7 @@ def _plot_sensors_3d( f"scales for {ch_type} must contain only numerical values, " f"got {scales} instead." ) - + ch_names = np.array(ch_names_all[ch_type], dtype="U") this_alpha = sensor_alpha[ch_type] if isinstance(sens_loc[0], dict): # meg coil if len(colors) == 1: @@ -1682,7 +1688,7 @@ def _plot_sensors_3d( template = sens_loc[idxs[0]] positions = np.array([sens_loc[i]["position"] for i in idxs]) quats = np.array([sens_loc[i]["quat"] for i in idxs]) - actor, _ = renderer.instanced_mesh( + actor, cloud = renderer.instanced_mesh( rr=template["rr"], tris=template["tris"], positions=positions, @@ -1691,6 +1697,7 @@ def _plot_sensors_3d( backface_culling=False, # visible from all sides ) actors[ch_type].append(actor) + cloud.field_data["ch_names"] = ch_names[idxs] else: # One GPU-instanced actor regardless of how many distinct # colors/scales are requested (broadcasting handles 1-vs-N). @@ -1720,7 +1727,7 @@ def _plot_sensors_3d( ) backface_culling = True actor_key = "eeg" - actor = _plot_glyphs( + actor, cloud = _plot_glyphs( renderer=renderer, loc=loc * unit_scalar, colors=these_colors, @@ -1737,6 +1744,7 @@ def _plot_sensors_3d( nearest=nearest, ) actors[actor_key].append(actor) + cloud.field_data["ch_names"] = ch_names[mask] actors = dict(actors) # get rid of defaultdict diff --git a/mne/viz/evoked_field.py b/mne/viz/evoked_field.py index e7eb16ce787..3f658485b8c 100644 --- a/mne/viz/evoked_field.py +++ b/mne/viz/evoked_field.py @@ -68,6 +68,10 @@ class EvokedField: The number of contours. .. versionadded:: 0.21 + contour_line_width : float + The line_width of the contour lines. + + .. versionadded:: 1.12 show_density : bool Whether to draw the field density as an overlay on top of the helmet/head surface. Defaults to ``True``. @@ -90,6 +94,16 @@ class EvokedField: ``True`` if there is more than one time point and ``False`` otherwise. .. versionadded:: 1.6 + background : tuple(int, int, int) + The color definition of the background: (red, green, blue). + + .. versionadded:: 1.12 + foreground : matplotlib color + Color of the foreground (will be used for colorbars and text). + None (default) will use black or white depending on the value + of ``background``. + + .. versionadded:: 1.12 %(verbose)s Notes @@ -112,11 +126,14 @@ def __init__( fig=None, vmax=None, n_contours=21, + contour_line_width=1, show_density=True, alpha=None, interpolation="nearest", interaction="terrain", time_viewer="auto", + background="black", + foreground=None, verbose=None, ): from .backends.renderer import _get_3d_backend, _get_renderer @@ -133,6 +150,7 @@ def __init__( self._vmax = _validate_type(vmax, (None, "numeric", dict), "vmax") self._n_contours = _ensure_int(n_contours, "n_contours") + self._contour_line_width = contour_line_width self._time_interpolation = _check_option( "interpolation", interpolation, @@ -141,6 +159,10 @@ def __init__( self._interaction = _check_option( "interaction", interaction, ["trackball", "terrain"] ) + self._bg_color = _to_rgb(background, name="background") + if foreground is None: + foreground = "w" if sum(self._bg_color) < 2 else "k" + self._fg_color = _to_rgb(foreground, name="foreground") surf_map_kinds = [surf_map["kind"] for surf_map in surf_maps] if vmax is None: @@ -192,13 +214,10 @@ def __init__( "is currently not supported inside a notebook." ) else: - self._renderer = _get_renderer( - fig, bgcolor=(0.0, 0.0, 0.0), size=(600, 600) - ) + self._renderer = _get_renderer(fig, bgcolor=background, size=(600, 600)) self._in_brain_figure = False self._units = "m" - self.plotter = self._renderer.plotter self.interaction = interaction # Prepare the surface maps @@ -236,7 +255,7 @@ def current_time_func(): if "%" in time_label: time_label = time_label % np.round(1e3 * time) self._time_label_actor = self._renderer.text2d( - x_window=0.01, y_window=0.01, text=time_label + x_window=0.01, y_window=0.01, text=time_label, color=foreground ) self._configure_dock() @@ -249,6 +268,10 @@ def current_time_func(): self._renderer.set_camera(azimuth=10, elevation=60, distance="auto") self._renderer.show() + @property + def plotter(self): + return self._renderer.plotter + def _prepare_surf_map(self, surf_map, color, alpha): """Compute all the data required to render a fieldlines map.""" from scipy.interpolate import interp1d @@ -325,6 +348,7 @@ def _prepare_surf_map(self, surf_map, color, alpha): vmin=-map_vmax, vmax=map_vmax, colormap=self._colormap_lines, + width=self._contour_line_width, ) else: contours = None # noqa @@ -370,7 +394,7 @@ def _update(self): if "%" in self._time_label: time_label = self._time_label % np.round(1e3 * self._current_time) self._time_label_actor = self._renderer.text2d( - x_window=0.01, y_window=0.01, text=time_label + x_window=0.01, y_window=0.01, text=time_label, color=self._fg_color ) self._renderer._update() @@ -416,6 +440,10 @@ def _callback(vmax, kind, scaling): ) r._layout_add_widget(layout, hlayout) + @_auto_weakref + def _rescale(): + self._rescale() + hlayout = r._dock_add_layout(vertical=False) r._dock_add_label( value="Rescale", @@ -424,19 +452,36 @@ def _callback(vmax, kind, scaling): ) r._dock_add_button( name="↺", - callback=self._rescale, + callback=_rescale, layout=hlayout, style="toolbutton", ) r._layout_add_widget(layout, hlayout) + @_auto_weakref + def _set_contours(n_contours): + self.set_contours(n_contours) + self._widgets["contours"] = r._dock_add_spin_box( name="Contour lines", value=21, rng=[0, 99], step=1, double=False, - callback=self.set_contours, + callback=_set_contours, + layout=layout, + ) + + @_auto_weakref + def _set_contour_line_width(line_width): + self.set_contour_line_width(line_width) + + self._widgets["contour_line_width"] = r._dock_add_slider( + name="Thickness", + value=self._contour_line_width, + rng=[0, 10], + callback=_set_contour_line_width, + double=True, layout=layout, ) r._dock_finalize() @@ -500,9 +545,13 @@ def _on_contours(self, event): break surf_map["contours"] = event.contours self._n_contours = len(event.contours) + if event.line_width is not None: + self._contour_line_width = event.line_width with disable_ui_events(self): if "contours" in self._widgets: self._widgets["contours"].set_value(len(event.contours)) + if "contour_line_width" in self._widgets and event.line_width is not None: + self._widgets["contour_line_width"].set_value(event.line_width) self._update() def set_time(self, time): @@ -537,6 +586,7 @@ def set_contours(self, n_contours): contours=np.linspace( -surf_map["map_vmax"], surf_map["map_vmax"], n_contours ).tolist(), + line_width=self._contour_line_width, ), ) @@ -571,3 +621,14 @@ def _rescale(self): current_data = surf_map["data_interp"](self._current_time) vmax = float(np.max(current_data)) self.set_vmax(vmax, kind=surf_map["map_kind"]) + + def set_contour_line_width(self, line_width): + """Set the line_width of the contour lines. + + Parameters + ---------- + line_width : float + The desired line_width of the contour lines. + """ + self._contour_line_width = line_width + self.set_contours(self._n_contours) diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 05694a109d5..d9b23fd6bcd 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -52,6 +52,7 @@ plot_sparse_source_estimates, set_3d_view, snapshot_brain_montage, + ui_events, ) from mne.viz._3d import _get_map_ticks, _linearize_map, _process_clim from mne.viz.utils import _fake_click, _fake_keypress, _fake_scroll, _get_cmap @@ -219,15 +220,27 @@ def test_plot_evoked_field(renderer): ) renderer.backend._close_all() - # Test some methods - fig = evoked.plot_field(maps, time_viewer=True) - assert isinstance(fig, EvokedField) + # Test some methods. Not all parameters are exposed through `plot_field`, so + # construct the `EvokedField` object directly. + fig = EvokedField( + evoked, + maps, + time_viewer=True, + contour_line_width=2, + background="white", + foreground="black", + ) + assert fig._contour_line_width == 2 + assert fig._widgets["contour_line_width"].get_value() == 2 fig._rescale() fig.set_time(0.05) assert fig._current_time == 0.05 fig.set_contours(10) assert fig._n_contours == 10 assert fig._widgets["contours"].get_value() == 10 + fig.set_contour_line_width(3) + assert fig._contour_line_width == 3 + assert fig._widgets["contour_line_width"].get_value() == 3 fig.set_vmax(2e-12, kind="meg") assert fig._surf_maps[1]["contours"][-1] == 2e-12 assert ( @@ -235,6 +248,18 @@ def test_plot_evoked_field(renderer): == DEFAULTS["scalings"]["grad"] * 2e-12 ) + # The contours (and their line width) can also be set through a UI event. + contours = [-2e-12, 0, 2e-12] + ui_events.publish( + fig, ui_events.Contours("field_strength_meg", contours, line_width=4) + ) + assert fig._n_contours == 3 + assert fig._contour_line_width == 4 + ui_events.publish( + fig, ui_events.Contours("field_strength_meg", contours, line_width=None) + ) + assert fig._contour_line_width == 4 # line_width=None keeps the current value + fig = evoked.plot_field(maps, time_viewer=False) assert isinstance(fig, Figure3D) renderer.backend._close_all() diff --git a/mne/viz/tests/test_ui_events.py b/mne/viz/tests/test_ui_events.py index 9fac7041999..3d64c8aae08 100644 --- a/mne/viz/tests/test_ui_events.py +++ b/mne/viz/tests/test_ui_events.py @@ -2,6 +2,8 @@ # License: BSD-3-Clause # Copyright the MNE-Python contributors. +from functools import partial + import matplotlib.pyplot as plt import pytest @@ -93,6 +95,23 @@ def callback(event): assert len(event_channels) == 0 +def test_subscriber_order(event_channels): + """Test that subscribers are called in the order in which they subscribed.""" + callback_calls = list() + + def callback(event, index): + callback_calls.append(index) + + fig = plt.figure() + callbacks = [partial(callback, index=index) for index in range(10)] + for cb in callbacks: + ui_events.subscribe(fig, "time_change", cb) + # Re-subscribing an existing callback keeps its original position. + ui_events.subscribe(fig, "time_change", callbacks[0]) + ui_events.publish(fig, ui_events.TimeChange(time=10.2)) + assert callback_calls == list(range(10)) + + def test_unsubscribe(event_channels): """Test unsubscribing from UI events.""" callback1_calls = list() diff --git a/mne/viz/ui_events.py b/mne/viz/ui_events.py index 825f301135b..25df95b6e2f 100644 --- a/mne/viz/ui_events.py +++ b/mne/viz/ui_events.py @@ -206,6 +206,9 @@ class Contours(UIEvent): kinds. contours : list of float The new values at which contour lines need to be drawn. + line_width : float | None + The line_width with which to draw the contour lines. Can be ``None`` to + indicate to keep using the current line_width. Attributes ---------- @@ -216,10 +219,14 @@ class Contours(UIEvent): kinds. contours : list of float The new values at which contour lines need to be drawn. + line_width : float | None + The line_width with which to draw the contour lines. Can be ``None`` to + indicate to keep using the current line_width. """ kind: str contours: list[str] + line_width: float | None = None @dataclass @@ -271,10 +278,10 @@ def _get_event_channel(fig): Returns ------- - channel : dict[event -> list] - The event channel. An event channel is a list mapping string event - names to a list of callback representing all subscribers to the - channel. + channel : dict[event -> dict] + The event channel. An event channel is a dict mapping string event + names to a dict of callbacks (used as an ordered set) representing all + subscribers to the channel, in the order in which they subscribed. """ import matplotlib @@ -342,7 +349,7 @@ def publish(fig, event, *, verbose=None): logger.debug(f"Publishing {event} on channel {fig}") for channel in channels: if event.name not in channel: - channel[event.name] = set() + channel[event.name] = dict() for callback in channel[event.name]: callback(event=event) @@ -360,12 +367,18 @@ def subscribe(fig, event_name, callback, *, verbose=None): callback : callable The function that should be called whenever the event is published. %(verbose)s + + Notes + ----- + Subscribers are called in the order in which they subscribed when the event + is published. """ channel = _get_event_channel(fig) logger.debug(f"Subscribing to channel {channel}") if event_name not in channel: - channel[event_name] = set() - channel[event_name].add(callback) + channel[event_name] = dict() + # use a dict as an ordered set: subscribers are called in subscription order + channel[event_name][callback] = None @verbose @@ -411,7 +424,7 @@ def unsubscribe(fig, event_names, callback=None, *, verbose=None): # Unsubscribe specific callback function. subscribers = channel[event_name] if callback in subscribers: - subscribers.remove(callback) + del subscribers[callback] else: warn( f'Cannot unsubscribe {callback} from event "{event_name}" ' diff --git a/tools/vulture_allowlist.py b/tools/vulture_allowlist.py index 31668930e5f..1fa94c4a34e 100644 --- a/tools/vulture_allowlist.py +++ b/tools/vulture_allowlist.py @@ -41,6 +41,7 @@ verbose_debug metadata_routing Plotter +customEvent # Decoding _.multi_class @@ -152,6 +153,8 @@ _._nearest_transformed_high_res_mri_idx_rpa _._nearest_transformed_high_res_mri_idx_nasion _._nearest_transformed_high_res_mri_idx_lpa +_.prop.culling +_.prop.lighting # Figures (prevent GC for example) _.decim_data diff --git a/tutorials/inverse/21_interactive_dipole_fit.py b/tutorials/inverse/21_interactive_dipole_fit.py new file mode 100644 index 00000000000..62ce94d449a --- /dev/null +++ b/tutorials/inverse/21_interactive_dipole_fit.py @@ -0,0 +1,185 @@ +""" +.. _tut-xfit: + +===================================================================== +Source localization by guided equivalent current dipole (ECD) fitting +===================================================================== + +This combination of manual specification and automated fitting is one of the oldest MEG +source estimation techniques :footcite:`Sarvas1987`. We will manually identify where and +when dipole source are active, upon which the fitting algorithm will find the best +location for the source. The result is a sparse source estimate of several equivalent +current dipoles (ECDs) that together explain (most of) the MEG evoked response. ECDs are +especially suited for capturing individual components of an evoked response (e.g. N100m, +N400m, etc.). Once the set of ECDs has been established, their timecourses can be +computed for multiple :class:`~mne.Evoked` objects, for example different experimental +conditions. + +This tutorial will demonstrate how to fit ECDs using the interactive GUI and also how to +drive that same GUI from Python code. Every screenshot below is of the *same* GUI +window, so you can follow along how its state evolves as we go. +""" + +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +# sphinx_gallery_preserve_gui = True + +# %% +# Guided ECD fitting using the GUI +# -------------------------------- +# +# Starting the GUI +# ~~~~~~~~~~~~~~~~ +# We can start the GUI from the command line by passing the ``mne dipolefit`` program +# an evoked file (filename typically ends in ``*-ave.fif``): +# +# .. code-block:: console +# +# $ mne dipolefit sample_audvis-ave.fif +# +# The GUI can also be started from an interactive python console. The minimal setup that +# can be used to fit dipoles is an :class:`~mne.Evoked` object and nothing else: +# +# .. code-block:: python +# +# mne.gui.dipolefit(evoked) +# +# When only specifying an evoked object, the GUI shows the sensors, a spherical head +# model and the electro-magnetic field recorded by the sensors, using an ad-hoc noise +# covariance matrix. If we provide more information, we can create a more accurate head +# model that provides better ECD fits and gives us more guidance for determining +# sources. On the command line there are various options you can use to specify files +# containing the covariance matrix, BEM model and MRI<->head transformation, see the +# output of ``mne dipolefit --help``. In an interactive python console, we can provide +# the appropriate MNE-Python objects when starting the GUI: + +import mne + +path = mne.datasets.sample.data_path() +meg_dir = path / "MEG" / "sample" +subjects_dir = path / "subjects" + +evoked = mne.read_evokeds(meg_dir / "sample_audvis-ave.fif", condition="Left Auditory") +evoked.apply_baseline() + +cov = mne.read_cov(meg_dir / "sample_audvis-cov.fif") +bem = mne.read_bem_solution( + subjects_dir / "sample" / "bem" / "sample-5120-5120-5120-bem-sol.fif" +) +trans = mne.read_trans(meg_dir / "sample_audvis_raw-trans.fif") + +# A distributed source estimate is a helpful guide for our dipole fits. +inv = mne.minimum_norm.read_inverse_operator( + meg_dir / "sample_audvis-meg-oct-6-meg-inv.fif" +) +stc = mne.minimum_norm.apply_inverse(evoked, inv) + +# Open the GUI with a better head model. +fitting_gui = mne.gui.dipolefit( + evoked, + cov=cov, + bem=bem, + trans=trans, + stc=stc, + ch_type="meg", # only use MEG sensors for this tutorial + subject="sample", + subjects_dir=subjects_dir, +) + +# %% +# Fitting a dipole +# ~~~~~~~~~~~~~~~~ +# During guided ECD fitting, we look for patterns in the electro-magnetic field to +# identify when and where sources may be active. We can use the time slider to examine +# how the field changes over time. The sample data is an evoked response to an auditory +# tone being played to the left of the participant and we can see the initial auditory +# response peaking at around 85 ms on the right hemisphere. The field shows a typical +# di-polar pattern with a pair of red/blue focii on either side of the source that +# should be located in auditory cortex (the distributed source estimate shows where it +# is). +# +# By pressing the "Fit dipole" button we instruct the algorithm to fit a dipole at the +# current time. After a few seconds of computation, the resulting dipole will be +# displayed as an arrow in the brain, indicating its source, as well as an arrow on the +# MEG helmet indicating the fit between the dipole and the field pattern. The timecourse +# of the dipole is shown below. On the right are controls to name, remove, temporarily +# (de-)activate, and save the dipole to a file. You also find a toggle switch to make +# the dipole's orientation dynamic or keep it fixed at the orientation it had at the +# time when it was fitted. +# +# Dragging the time slider and clicking the "Fit dipole" button have programmatic +# equivalents, which is what we will use here to build up the figures in this tutorial: + +fitting_gui.set_time(0.085) +fitting_gui.fit_dipole() + +# %% +# Selecting channels to guide the ECD modeling +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# At nearly the same time, auditory responses are occurring in both left and right +# auditory cortex. Hence, the single dipole that we fitted without any guidance will +# have some bias as the algorithm attempted to fit the ECD to the entire bi-lateral +# field pattern. +# +# To isolate portions of the field pattern that contain a single pair of red/blue focii, +# the ideal fitting target for the algorithm, we can restrict the analysis to a subset +# of sensors. To do so, first press the "Sensor data" button, which will open a new +# window showing the evoked response across all sensors. By clicking and dragging the +# mouse we can make a lasso selection around the sensors we wish to include in the +# analysis. Hold ``CTRL`` to add to the current selection and ``CTRL + SHIFT`` to remove +# from the current selection. The currently selected sensors are highlighted in green in +# the main window, showing the portion of the field pattern they cover. When you are +# happy with the selection, you can use the "Fit dipole" button as before to fit a +# dipole using the selected sensors at the current timepoint. +# +# Remove or de-activate the dipole we previously fitted to the entire field pattern and +# fit two dipoles using the left-side and right-side sensors respectively. It is helpful +# to name them. This part of the workflow is inherently interactive, so we cannot +# reproduce it here in code. + +# %% +# Multi/single dipole modes +# ~~~~~~~~~~~~~~~~~~~~~~~~~ +# By default, the fitting algorithm is in "Multi dipole (MNE)" mode, meaning portions of +# the signal attributed to one dipole can not be attributed to a second dipole at the +# same time. You will notice that if you have two dipoles with similar orientations +# close to each other, their timecourses become a strange mixture as each dipole will +# claim a part of the same signal. To prevent this, we can switch the algorithm over to +# "Single dipole" using the "Dipole model" dropdown. In this mode, the timecourse of +# each dipole will be computed whilst ignoring all other dipoles, which is useful when +# evaluating multiple candidate dipoles for the same source. + +# %% +# Saving and loading sets of dipoles +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# We can save the fitted dipoles using the "Save dipoles" button. There are two possible +# file formats for this: plain text (``.dip``) and a binary format (``.bdip``). These +# formats are compatible with MEGIN's software, allowing interoperability between +# MNE-Python and Xfit. Saved dipoles can be read back with :func:`mne.read_dipole` and +# added to an existing dipole fitting GUI, optionally under a name of our choosing. +# Here, we add the dipole that a lasso selection of the left-side sensors would have +# produced: + +dips_to_add = mne.read_dipole(meg_dir / "sample_audvis_set1.dip") +dips_to_add = dips_to_add[[33]] # add only one of the 34 dipoles in the file +fitting_gui.add_dipole(dips_to_add, name="lh") + +# %% +# Working with the dipoles from Python +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Although each dipole is fitted at a single time point, its timecourse is estimated +# across the entire epoch, so we can move the time cursor to any latency to inspect the +# dipole model there. The dipoles themselves are available as a list of +# :class:`mne.Dipole` objects and can be saved without touching the GUI at all: + +fitting_gui.set_time(0.115) +fitted_dipoles = fitting_gui.dipoles # the dipoles we fitted +print([dip.name for dip in fitted_dipoles]) +# save with: fitting_gui.save("my_file.dip") + +# %% +# References +# ---------- +# .. footbibliography::