restrict_forward_to_label vs stc.in_label methods

[dasdiptyajit]

Description of the problem

I am testing sensitivity profiles/maps of different sensor types ('grad', 'eeg', 'mag') for different labels.

During testing, I tried 2 methods:

Method1: restrict the forward sol' with some particular labels and then compute sensitivity maps (mean across vertices of the labels).

Method2: compute sensitivity maps for a full forward sol' and use stc.in_label () method by passing the labels to extract the particular labels sensitivity.

Expected outcomes: since both of the methods follow similar ways of calculating sensitivity maps. The outcomes of these two should be similar.

Steps to reproduce

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
=================================================================
test restrict_forward_to_label vs stc.in_label method
on sensitivity profile using labels
=================================================================
"""

print(__doc__)

# --------------------------------------
# import modules
# --------------------------------------
from os.path import join
import mne
from mne.datasets import sample
from mne import make_forward_solution

# load data
data_path = sample.data_path()
subject = 'sample'
subjects_dir = join(data_path, 'subjects')
meg_path = join(data_path, 'MEG', 'sample')

# read evoked
evoked_fname = join(meg_path, 'sample_audvis-ave.fif')
evoked = mne.read_evokeds(evoked_fname)[0]
evoked.pick_types(meg=True, eeg=True).apply_baseline((None, 0))

# create a forward sol
bem = join(data_path,  'subjects', 'sample', 'bem', 'sample-5120-5120-5120-bem-sol.fif')
src = join(data_path,  'subjects', 'sample', 'bem', 'sample-oct-6-src.fif')
trans = join(meg_path, 'sample_audvis_raw-trans.fif')

# make forward solution
fwd = make_forward_solution(evoked.info, trans, src, bem, meg=True, eeg=True, mindist=3.0, n_jobs=3)

# convert the forward solution with surface orientation /normal orientation
fwd = mne.convert_forward_solution(fwd, surf_ori=True, copy=False, force_fixed=True)

# method 1: restrict_forward_to_label
# read labels
labels = mne.read_labels_from_annot("fsaverage", "HCPMMP1", subjects_dir=subjects_dir)
aud_label_lh = [label for label in labels if label.name == "L_A1_ROI-lh"][0]
aud_label_rh = [label for label in labels if label.name == "R_A1_ROI-rh"][0]

# morphed the labels for individual subject anatomy
subject_labels = mne.morph_labels([aud_label_lh, aud_label_rh], subject_to=subject, subject_from='fsaverage',
                                  subjects_dir=subjects_dir, verbose=None)
# restrict fwd to auditory labels
fwd_aud_labels = mne.forward.restrict_forward_to_label(fwd.copy(), subject_labels)

# test on different sensors
ch_types = ['grad', 'eeg', 'mag']

for ch_type in ch_types:
    # compute sensitivity maps with forward labels
    grad_map_m1 = mne.sensitivity_map(fwd_aud_labels, ch_type=ch_type, mode='fixed')
    m1_mean = grad_map_m1.data.mean()  # mean sensitivity of the labels

    # method 2: 'in_label' method

    # compute sensitivity maps with full forward
    grad_map_full = mne.sensitivity_map(fwd, ch_type=ch_type, mode='fixed')

    # now compute the sensitivity for auditory labels using stc: 'in_label' method
    grad_map_m2 = grad_map_full.in_label(subject_labels[0]+subject_labels[1])

    m2_mean = grad_map_m2.data.mean()  # mean sensitivity of the labels

    print('%s method1 output:' % ch_type, round(m1_mean, 2))
    print('%s method2 output:' % ch_type, round(m2_mean, 2))

Link to data

No response

Expected results

The outcomes of these two methods should be similar.

Actual results

grad method1 output: 0.69
grad method2 output: 0.26

eeg method1 output: 0.86
eeg method2 output: 0.46

mag method1 output: 0.74
mag method2 output: 0.36

Additional information

mne sys_info() output

Platform Linux-5.19.0-41-generic-x86_64-with-glibc2.35
Python 3.9.16 (main, Mar 1 2023, 18:22:10) [GCC 11.2.0]
Executable /home/dip_meg/anaconda3/envs/mne/bin/python
CPU x86_64 (4 cores)
Memory 15.5 GB

Core
├☑ mne 1.4.0
.
.

[dasdiptyajit]

@larsoner can you have a look at this.

thanks,
Dip

[agramfort]

no they should not be the same as the inverse operator sees "all available sources" when doing the unmixing. Can this explain the issue?

[dasdiptyajit]

@agramfort Hi!

Isn't it the forward matrix that we are only considering here to calculate the sensitivity maps? I don't see any unmixing/inverse in the original mne code. As per I know, the sensitivity maps/topography for each dipole is the column of the leadfield matrix. I would expect the sensitivity topography remains intact over a type sensors for individual dipole. It should be independent of computing a full brain or small ROI forward sol'.

Or, are you saying computing a forward| Ledfields for a dipole (say : source space with single dipole) will be different than for a full brain forward| Leadfields and extracting only the forward| Leadfields for that dipole? In that case, creating these maps with a small ROI using restrict_forward_to_label would be overestimated?

thanks!

[larsoner]

I think it's due to the normalization by the max:

mne-python/mne/proj.py

Lines 512 to 513 in 68e4f89

	if mode in ["fixed", "free"]:
	sensitivity_map /= np.max(sensitivity_map)

We should probably add a note that for 'fixed' and 'free' methods this normalization happens

[larsoner]

Argh no there is some deeper bug. I'll keep looking...

[dasdiptyajit]

Yes! Even without the normalization, values are different between M1 and M2.