Noisy resolution matrix

doc/changes/dev/11092.newfeature.rst

@@ -0,0 +1 @@
+- Added noisy point-spread functions to mne.minimum_norm.make_inverse_resolution_matrix (:gh:`11058` by `Olaf Hauk`_)

examples/inverse/psf_empirical.py

@@ -0,0 +1,147 @@
+"""
+.. _ex-psf-empirical:
+
+===================================================
+Plot point-spread functions (PSFs) with added noise
+===================================================
+
+Visualise the point-spread of a vertex in a sLORETA source estimate when we assume the
+data to be noisy.
+"""
+# Authors: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
+#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
+#
+# License: BSD-3-Clause
+# Copyright the MNE-Python contributors.
+
+import mne
+from mne.datasets import sample
+from mne.minimum_norm import get_point_spread, make_inverse_resolution_matrix
+
+print(__doc__)
+
+data_path = sample.data_path()
+subjects_dir = data_path / "subjects"
+meg_path = data_path / "MEG" / "sample"
+fname_fwd = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
+fname_cov = meg_path / "sample_audvis-cov.fif"
+fname_evo = meg_path / "sample_audvis-ave.fif"
+
+# read forward solution
+forward = mne.read_forward_solution(fname_fwd)
+# forward operator with fixed source orientations
+mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True, copy=False)
+
+# noise covariance matrix
+noise_cov = mne.read_cov(fname_cov)
+
+# evoked data for info
+info = mne.io.read_info(fname_evo)
+
+# make inverse operator from forward solution
+# free source orientation
+inverse_operator = mne.minimum_norm.make_inverse_operator(
+    info=info, forward=forward, noise_cov=noise_cov, loose=0.0, depth=None
+)
+
+# regularisation parameter
+snr = 3.0
+lambda2 = 1.0 / snr**2
+method = "MNE"
+
+# %%
+# We compute noisy PSFs for two locations, one deep and one superficial.
+# For each location we compute two PSFs::
+# 1) For a high SNR, approximating the noise-less case.
+# 2) With the SNR of regularization and the added noise matched (realistic).
+
+# vertices of a deep and a superficial source
+sources = {"deep": 2146, "super": 2209}
+
+# point-spread functions of the deep and superficial source
+stc_psf = {"deep": {}, "super": {}}
+
+# for a deep source
+# compute resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["deep"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["deep"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["deep"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["deep"]], norm=True
+)
+del rm_highsnr
+
+# the same for a superficial source
+# compute matching resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["super"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["super"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["super"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["super"]], norm=True
+)
+del rm_highsnr
+
+##############################################################################
+# Visualize
+# ---------
+# PSF:
+
+vertno_lh = forward["src"][0]["vertno"]
+for source in sources:
+    # Which vertex corresponds to selected source
+    verttrue = [vertno_lh[sources[source]]]  # just one vertex
+    for psf in ["matched", "highsnr"]:
+        brain_psf = stc_psf[source][psf].plot(
+            "sample", "inflated", "lh", subjects_dir=subjects_dir
+        )
+        brain_psf.show_view("ventral")
+        title_str = f"{method} {source} {psf}"
+        brain_psf.add_text(
+            0.1, 0.9, title_str, "title", font_size=16, color=(255, 255, 255)
+        )
+        # mark true source location
+        brain_psf.add_foci(
+            verttrue, coords_as_verts=True, scale_factor=1.0, hemi="lh", color="green"
+        )
+
+# %%
+# The green spheres indicate the true source location.
+# For high SNR, PSFs for neither location are ideal. For the superficial
+# it is at least centred around the true source location, but still widespread
+# and with sidelobes. For the deep location the peaks occur at large distances
+# to the true location. This shows that spatial resolution is not an issue of
+# insufficient SNR but is inherently limited.
+# For matched (more realistic) SNR the PSF resembles the noise-less PSF, albeit
+# with more widely spread-out sidelobes. This indicates that moderate
+# regularization still achieves good spatial resolution while effectively
+# suppressing noise.

mne/minimum_norm/resolution_matrix.py

@@ -11,7 +11,7 @@
 from mne.minimum_norm.inverse import InverseOperator
 
 from .._fiff.constants import FIFF
-from .._fiff.pick import pick_channels_forward
+from .._fiff.pick import pick_channels_forward, pick_types
 from ..evoked import EvokedArray
 from ..forward.forward import Forward, convert_forward_solution
 from ..label import Label
@@ -24,12 +24,19 @@
     logger,
     verbose_static,
 )
-from .inverse import apply_inverse
+from .inverse import apply_inverse, apply_inverse_cov
+from .spatial_resolution import _rectify_resolution_matrix
 
 
 @verbose_static()
 def make_inverse_resolution_matrix(
-    forward, inverse_operator, method="dSPM", lambda2=1.0 / 9.0, verbose=None
+    forward,
+    inverse_operator,
+    method="dSPM",
+    lambda2=1.0 / 9.0,
+    noise_cov=None,
+    snr=None,
+    verbose=None,
 ):
     """Compute resolution matrix for linear inverse operator.
 
@@ -43,6 +50,11 @@ def make_inverse_resolution_matrix(
         Inverse method to use (MNE, dSPM, sLORETA).
     lambda2 : float
         The regularisation parameter.
+    noise_cov : None | instance of Covariance
+        Noise covariance matrix to compute noise power in source space.
+    snr : None | float
+        Signal-to-noise ratio for source signal vs noise power.
+        If None, SNR is inferred from lambda2 (as 1/sqrt(lambda2)).
     verbose : bool | str | int | None
         Control verbosity of the logging output. If ``None``, use the default
         verbosity level. See the :ref:`logging documentation <tut-logging>` and
@@ -52,13 +64,29 @@ def make_inverse_resolution_matrix(
     Returns
     -------
     resmat: array, shape (n_orient_inv * n_dipoles, n_orient_fwd * n_dipoles)
-        Resolution matrix (inverse operator times forward operator).
-        The result of applying the inverse operator to the forward operator.
+        Resolution matrix (inverse operator times forward operator). The
+        columns of the resolution matrix are the point-spread functions (PSFs)
+        and the rows are the cross-talk functions (CTFs).
         If source orientations are not fixed, all source components will be
         computed (i.e. for n_orient_inv > 1 or n_orient_fwd > 1).
-        The columns of the resolution matrix are the point-spread functions
-        (PSFs) and the rows are the cross-talk functions (CTFs).
+        If 'noise_cov' and 'snr' are not None, then source noise estimated from
+        the noise covariance matrix will be added to columns of the resolution
+        matrix (PSFs). In this case, 'resmat' will be of shape
+        (n_dipoles, n_orient_fwd * n_dipoles) with values pooled across
+        n_orient_inv source orientations per location.
+        Note: If the resolution matrix is computed with a noise covariance
+        matrix then only its columns, i.e. PSFs, can meaningfully be
+        interpreted. It must not be used to compute CTFs or resolution metrics
+        for CTFs.
     """
+    if noise_cov is None and snr is not None:
+        msg = "snr should be None if noise_cov is None."
+        raise ValueError(msg)
+
+    if noise_cov is not None and snr is None:
+        snr = np.sqrt(1.0 / lambda2)
+        logger.info(f"Inferring snr from lamdba2 as: {snr}.")
+
     # make sure forward and inverse operator match
     inv = inverse_operator
     fwd = _convert_forward_match_inv(forward, inv)
@@ -77,6 +105,55 @@ def make_inverse_resolution_matrix(
     logger.info(
         f"Dimensions of resolution matrix: {resmat.shape[0]} by {resmat.shape[1]}."
     )
+
+    # add source noise power to columns of resolution matrix
+    if noise_cov is not None:
+        info = _prepare_info(inv)
+        # compute source noise power
+        stc = apply_inverse_cov(
+            noise_cov,
+            info,
+            inv,
+            nave=1,
+            lambda2=lambda2,
+            method=method,
+            pick_ori=None,
+            prepared=False,
+            label=None,
+            method_params=None,
+            use_cps=True,
+            verbose=None,
+        )
+
+        # output will be square-root of power, so take intensities before
+        # adding noise
+        shape = resmat.shape
+        if not shape[0] == shape[1]:
+            # pool loose/free orientations
+            resmat = _rectify_resolution_matrix(resmat)
+        else:
+            resmat = np.abs(resmat)
+
+        # scaling factor
+
+        idx_ch_types = {}
+        idx_ch_types["eeg"] = pick_types(info, eeg=True, meg=False)
+        idx_ch_types["mag"] = pick_types(info, eeg=False, meg="mag")
+        idx_ch_types["gra"] = pick_types(info, eeg=False, meg="grad")
+
+        # scaling similar to Samuelsson et al., Neuroimage 2021 (p. 4)
+        alphas = {}
+        for cht in idx_ch_types:
+            # signal intensity of all sources across channels of this type
+            alphas[cht] = np.sqrt((leadfield[idx_ch_types[cht], :] ** 2).mean())
+            # divided by noise intensity across channels of this type
+            alphas[cht] /= np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean())
+
+        alpha = (3 * snr) / (alphas["eeg"] + alphas["mag"] + alphas["gra"])
+
+        # Add square root of source noise power to every column, scale
+        # depending on SNR
+        resmat = resmat + (1 / alpha) * np.sqrt(stc.data)
     return resmat
 
 
@@ -94,7 +171,29 @@ def _get_psf_ctf(
     vector=False,
     verbose=None,
 ):
-    """Get point-spread (PSFs) or cross-talk (CTFs) functions."""
+    """Get point-spread (PSFs) or cross-talk (CTFs) functions.
+
+    Parameters
+    ----------
+    resmat : array, shape (n_orient_inv * n_dipoles, n_orient_fwd * n_dipoles)
+        Forward Operator.
+    src : Source Space
+        Source space used to compute resolution matrix.
+    %(idx_pctf)s
+    func : str ('psf' | 'ctf')
+        Whether to produce PSFs or CTFs. Defaults to psf.
+    %(mode_pctf)s
+    %(n_comp_pctf_n)s
+    %(norm_pctf)s
+    %(return_pca_vars_pctf)s
+    %(vector_pctf)s
+    %(verbose)s
+
+    Returns
+    -------
+    %(stcs_pctf)s
+    %(pca_vars_pctf)s
+    """
     # check for consistencies in input parameters
     _check_get_psf_ctf_params(mode, n_comp, return_pca_vars)
 
@@ -128,7 +227,7 @@ def _get_psf_ctf(
         (n_verts != n_c) and (n_c / 3 != n_verts)
     ):
         msg = (
-            f"Number of vertices ({n_verts}) and corresponding dimension of"
+            f"Number of vertices ({n_verts}) and corresponding dimension of "
             f"resolution matrix ({n_r}, {n_c}) do not match"
         )
         raise ValueError(msg)
@@ -181,7 +280,6 @@ def _get_psf_ctf(
                     funcs_vert = funcs[3 * i : 3 * i + 3, :]
                     funcs_int[i, :] = np.sqrt((funcs_vert**2).sum(axis=0))
                 funcs = funcs_int
-
         stc = _make_stc(
             funcs,
             vertno,
@@ -192,6 +290,7 @@ def _get_psf_ctf(
             vector=vector,
             source_nn=nn,
         )
+
         stcs.append(stc)
         pca_vars.append(pca_var)
 
@@ -223,7 +322,7 @@ def _vertices_for_get_psf_ctf(idx, src):
     """Get vertices in source space for PSFs/CTFs in _get_psf_ctf()."""
     # idx must be list
     # if label(s) specified get the indices, otherwise just carry on
-    if type(idx[0]) is Label:
+    if isinstance(idx[0], Label):
         # specify without source time courses, gets indices per label
         verts_labs, _ = _prepare_label_extraction(
             stc=None,
@@ -269,7 +368,7 @@ def _normalise_psf_ctf(funcs, norm):
     return funcs
 
 
-def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars, nn):
+def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
     """Summarise PSFs/CTFs across vertices."""
     s_var = None  # only computed for return_pca_vars=True
 
@@ -603,7 +702,7 @@ def _get_matrix_from_inverse_operator(
         The inverse operator.
     forward : instance of Forward
         The forward operator.
-    method : 'MNE' | 'dSPM' | 'sLORETA'
+    method : 'MNE' | 'dSPM' | 'sLORETA' | 'eLORETA'
         Inverse methods (for apply_inverse).
     lambda2 : float
         The regularization parameter (for apply_inverse).