MRG, ENH: Add forward sensitivity maps to reports

doc/changes/dev/14107.newfeature.rst

@@ -0,0 +1 @@
+Added the ``sensitivity`` parameter to :meth:`mne.Report.add_forward` to include forward sensitivity maps in reports, by :newcontrib:`Mariam Husain`.

doc/changes/names.inc

@@ -267,6 +267,7 @@
 .. _Manorama Kadwani: https://github.com/synapse-stream
 .. _Manu Sutela: https://github.com/MJAS1
 .. _Marcin Koculak: https://github.com/mkoculak
+.. _Mariam Husain: https://github.com/mariam-hedgie
 .. _Marian Dovgialo: https://github.com/mdovgialo
 .. _Marijn van Vliet: https://github.com/wmvanvliet
 .. _Mark Alexander Henney: https://github.com/henneysq

mne/report/report.py

@@ -41,8 +41,8 @@
 from ..minimum_norm import InverseOperator, read_inverse_operator
 from ..parallel import parallel_func
 from ..preprocessing.ica import read_ica
-from ..proj import read_proj
-from ..source_estimate import SourceEstimate, read_source_estimate
+from ..proj import read_proj, sensitivity_map
+from ..source_estimate import _BaseSourceEstimate, read_source_estimate
 from ..source_space._source_space import _ensure_src
 from ..surface import dig_mri_distances
 from ..transforms import _find_trans
@@ -89,6 +89,13 @@
 
 _BEM_VIEWS = ("axial", "sagittal", "coronal")
 
+# constant for mapping of internal channel type keys to better labels
+# used in figures
+_SENSITIVITY_MAP_CH_LABELS = (
+    ("grad", "Gradiometers"),
+    ("mag", "Magnetometers"),
+    ("eeg", "EEG"),
+)
 
 # For raw files, we want to support different suffixes + extensions for all
 # supported file formats
@@ -1583,6 +1590,7 @@ def add_forward(
         subject=None,
         subjects_dir=None,
         plot=False,
+        sensitivity=False,
         tags=("forward-solution",),
         section=None,
         replace=False,
@@ -1604,6 +1612,10 @@ def add_forward(
             If True, plot the source space of the forward solution.
 
             .. versionadded:: 1.10
+        sensitivity : bool
+            If True, render sensitivity maps for all available sensor types.
+
+            .. versionadded:: 1.13
         %(tags_report)s
         %(section_report)s
 
@@ -1626,6 +1638,7 @@ def add_forward(
             tags=tags,
             replace=replace,
             plot=plot,
+            sensitivity=sensitivity,
         )
 
     @fill_doc
@@ -3709,6 +3722,7 @@ def _add_forward(
         title,
         image_format,
         plot,
+        sensitivity,
         section,
         tags,
         replace,
@@ -3720,9 +3734,17 @@ def _add_forward(
         subject = self.subject if subject is None else subject
         subject = forward["src"][0]["subject_his_id"] if subject is None else subject
 
-        # XXX Todo
-        # Render sensitivity maps
         sensitivity_maps_html = ""
+        if sensitivity:
+            sensitivity_maps_html = self._render_forward_sensitivity_maps(
+                forward=forward,
+                subject=subject,
+                subjects_dir=subjects_dir,
+                image_format=image_format,
+                section=section,
+                title=title,
+                tags=tags,
+            )
         source_space_html = ""
         if plot:
             source_space_html = self._src_html(
@@ -3749,6 +3771,63 @@ def _add_forward(
             replace=replace,
         )
 
+    def _render_forward_sensitivity_maps(
+        self, *, forward, subject, subjects_dir, image_format, section, title, tags
+    ):
+        # render sensitivity maps for all available sensors as one slider
+        subjects_dir = self.subjects_dir if subjects_dir is None else subjects_dir
+
+        if get_3d_backend() is None:
+            raise RuntimeError(
+                "A 3D backend is needed to render source estimates in a report."
+            )
+        stc_plot_kwargs = _handle_default("report_stc_plot_kwargs", None)
+        stc_plot_kwargs.update(
+            subject=subject,
+            subjects_dir=subjects_dir,
+            clim=dict(kind="value", lims=(0, 0.5, 1.0)),
+        )
+        if self.img_max_width is not None:
+            stc_plot_kwargs["size"] = (
+                stc_plot_kwargs["size"][0],
+                min(stc_plot_kwargs["size"][1], self.img_max_width),
+            )
+        ch_types = forward["info"].get_channel_types(unique=True)
+        figs = []
+        captions = []
+        # compute each sensitivity map and plot it
+        for ch_type, label in _SENSITIVITY_MAP_CH_LABELS:
+            if ch_type not in ch_types:
+                continue
+            stc = sensitivity_map(forward, ch_type=ch_type)
+            brain = stc.plot(**stc_plot_kwargs)
+            brain._renderer.plotter.subplot(0, 0)
+            brain.add_text(
+                x=0.5,
+                y=0.8,
+                text=label,
+                font_size=8,
+                justification="center",
+            )
+            figs.append(brain.screenshot(time_viewer=True, mode="rgb"))
+            brain.close()
+            captions.append(label)
+            # make one working slider
+        html_partial = self._render_slider(
+            figs=figs,
+            imgs=None,
+            captions=captions,
+            title="Sensitivity maps",
+            start_idx=0,
+            image_format=image_format,
+            tags=tags,
+            klass="stc",
+            own_figure=False,
+        )
+        return html_partial(
+            id_=self._get_dom_id(section=section, title=f"{title}-sensitivity")
+        )
+
     def _src_html(
         self,
         *,
@@ -4508,7 +4587,38 @@ def _add_stc(
         replace,
     ):
         """Render STC."""
-        if isinstance(stc, SourceEstimate):
+        html_partial = self._render_stc(
+            stc=stc,
+            title=title,
+            subject=subject,
+            subjects_dir=subjects_dir,
+            n_time_points=n_time_points,
+            image_format=image_format,
+            tags=tags,
+            stc_plot_kwargs=stc_plot_kwargs,
+        )
+        self._add_or_replace(
+            title=title,
+            section=section,
+            tags=tags,
+            html_partial=html_partial,
+            replace=replace,
+        )
+
+    def _render_stc(
+        self,
+        *,
+        stc,
+        title,
+        subject,
+        subjects_dir,
+        n_time_points,
+        image_format,
+        tags,
+        stc_plot_kwargs,
+    ):
+        """Render an STC as embeddable report HTML."""
+        if isinstance(stc, _BaseSourceEstimate):
             if subject is None:
                 subject = self.subject  # supplied during Report init
                 if not subject:
@@ -4543,96 +4653,40 @@ def _add_stc(
             )
         t_zero_idx = np.abs(times).argmin()  # index of time closest to zero
 
-        # Plot using 3d backend if available, and use Matplotlib
-        # otherwise.
-        import matplotlib.pyplot as plt
+        if get_3d_backend() is None:
+            raise RuntimeError(
+                "A 3D backend is required to render source estimates in a report."
+            )
 
         stc_plot_kwargs = _handle_default("report_stc_plot_kwargs", stc_plot_kwargs)
         stc_plot_kwargs.update(subject=subject, subjects_dir=subjects_dir)
-        # we need to set the size based on the min (img_max_width can be None)
         if self.img_max_width is not None:
             stc_plot_kwargs["size"] = (
                 stc_plot_kwargs["size"][0],
                 min(stc_plot_kwargs["size"][1], self.img_max_width),
             )
-        if get_3d_backend() is not None:
-            brain = stc.plot(**stc_plot_kwargs)
-            brain._renderer.plotter.subplot(0, 0)
-            backend_is_3d = True
-        else:
-            backend_is_3d = False
+
+        brain = stc.plot(**stc_plot_kwargs)
+        brain._renderer.plotter.subplot(0, 0)
 
         figs = []
         for t in times:
-            with warnings.catch_warnings():
-                warnings.filterwarnings(
-                    action="ignore",
-                    message="More than 20 figures have been opened",
-                    category=RuntimeWarning,
-                )
-
-                if backend_is_3d:
-                    brain.set_time(t)
-                    figs.append(brain.screenshot(time_viewer=True, mode="rgb"))
-                else:
-                    fig_lh = plt.figure(layout="constrained")
-                    fig_rh = plt.figure(layout="constrained")
-
-                    brain_lh = stc.plot(
-                        views="lat",
-                        hemi="lh",
-                        initial_time=t,
-                        backend="matplotlib",
-                        subject=subject,
-                        subjects_dir=subjects_dir,
-                        figure=fig_lh,
-                    )
-                    brain_rh = stc.plot(
-                        views="lat",
-                        hemi="rh",
-                        initial_time=t,
-                        subject=subject,
-                        subjects_dir=subjects_dir,
-                        backend="matplotlib",
-                        figure=fig_rh,
-                    )
-                    _constrain_fig_resolution(
-                        fig_lh,
-                        max_width=stc_plot_kwargs["size"][0],
-                        max_res=self.img_max_res,
-                    )
-                    _constrain_fig_resolution(
-                        fig_rh,
-                        max_width=stc_plot_kwargs["size"][0],
-                        max_res=self.img_max_res,
-                    )
-                    figs.append(brain_lh)
-                    figs.append(brain_rh)
-                    plt.close(fig_lh)
-                    plt.close(fig_rh)
-
-        if backend_is_3d:
-            brain.close()
-        else:
-            brain_lh.close()
-            brain_rh.close()
+            brain.set_time(t)
+            figs.append(brain.screenshot(time_viewer=True, mode="rgb"))
+        brain.close()
 
         captions = [f"Time point: {round(t, 3):0.3f} s" for t in times]
-        self._add_slider(
+        return self._render_slider(
             figs=figs,
             imgs=None,
             captions=captions,
             title=title,
             image_format=image_format,
             start_idx=t_zero_idx,
-            section=section,
             tags=tags,
-            replace=replace,
-            own_figure=False,  # prevent rescaling
+            klass="stc",
+            own_figure=False,
         )
-        for fig in figs:
-            if not isinstance(fig, np.ndarray):
-                plt.close(fig)
 
     @_use_agg
     def _add_bem(