restrict_forward_to_label gives unusual 'vertno'

[clamus]

Hi Guys,

I noticed that restrict_forward_to_label returns a fwd where vertno = fwd['src'][0 of 1]['vertno'] are indices that do not correspond to the values of vertices in the original source space (like from 0 to 200k). Instead, vertno gives the vertices used for selecting the columns in the original forward solution data directly. Is this the desired behavior? Or should vertno have the indices relating the restricted forward solution to the source space surface?

[clamus]

I put together a tiny example to illustrate the issue---if this is an issue. This example takes evoked, covariance, forward data from mne's test data sets. Then restricts the forward to a label in the left auditory cortex. Computes the minimum norm inverse solution using the restricted forward, as well as the respective stc object for visualization. What one finds is that the solution is plotted on the visual cortex, which does not make sense since the forward model is restricted to a tiny label in left auditory cortex. The reason this happens is because restrict_forward_to_label creates a source space in which the field vertno in each hemisphere has index values that correspond to the columns of the original forward solution, and not to vertices in the source space surface. The question is then, if this is by design, or if it is an issue. If it is an issue, I can fix it :)

import os.path as op

from mne import read_evokeds, read_forward_solution, read_cov, read_label
from mne.datasets import testing
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne.forward import restrict_forward_to_label

# Get data paths
data_path = testing.data_path(download=False)
subjects_dir = op.join(data_path, 'subjects')
fwd_fname = op.join(data_path, 'MEG', 'sample',
                    'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif')
evoked_fname = op.join(data_path, 'MEG', 'sample',
                       'sample_audvis_trunc-ave.fif')
cov_fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-cov.fif')

# Read data
evoked = read_evokeds(evoked_fname, condition='Left Auditory',
                      baseline=(None, 0), proj=True)
fwd = read_forward_solution(fwd_fname, force_fixed=True)
cov = read_cov(cov_fname)
info = evoked.info

# Read auditory label in lh
label_aud_lh_fname = op.join(data_path, 'MEG', 'sample', 'labels',
                             'Aud-lh.label')
label_aud_lh = read_label(label_aud_lh_fname, subject='sample')

# Restrict fwd to auditory label
fwd_label = restrict_forward_to_label(fwd, label_aud_lh)

# Compute inverse on label
inv = make_inverse_operator(info, fwd_label, cov, fixed=True, depth=None)
method = "dSPM"
snr = 3.
lambda2 = 1. / snr ** 2
stc = apply_inverse(evoked, inv, lambda2, method=method, pick_ori=None)

# Plot
# NOTE: The estimates fall in visual cortex, while the label restricting
# the forward is in lh auditory cortex

brain = stc.plot(subject='sample', subjects_dir=subjects_dir, colormap='hot')
brain.scale_data_colormap(fmin=0, fmid=0.5, fmax=3, transparent=True)
brain.add_label(label_aud_lh, color='magenta')

[dgwakeman]

Looks like a bug to me, but I honestly don't understand the point of this function. I get that it does something, but this feels like feature creep to me. It can also be dangerous if people don't understand the implications of using it. In your example, it is basically just a dipole. this seems like something a user could accidentally use, and not realize how dangerous it is. my 2cents

[clamus]

Thanks for the input @dgwakeman. The reason I wanted to use restrict_forward_to_label was to make a kind of complex simulation of distributed sources with some particular dynamics over extended areas of cortex, but not all of cortex. In my case, the functionality of simulation.sources does not help me. What I do is to select some particular areas (labels), then simulate source dynamics outside mne, and then generate the MEG/EEG data using the forward model restricted to the labels. I am sure there are other ways to go about doing what I need as infinitely many road lead to Rome.

The question then is if this issue need to be fixed.

[dgwakeman]

I'll leave that to higher up devs than me, like @Eric89GXL, @dengemann, and @agramfort. FWIW, why doesn't simulation.sources help you? It seems pretty simple to just generate an STC of whatever you want, where you want it and then run that function on it. d

[clamus]

I went ahead and fixed this issue. I suspect the silence of the higher ups relates to the March 5 deadline for submitting many NIH grants.

[dgwakeman]

@clamus any reason simulation.sources doesn't work?

[larsoner]

I can't remember what restrict_forward_to_label was originally designed for, but it sounds like there is indeed a bug with the vertno. It would be nice if we could enhance the simulation module to do what you need, if as @dgwakeman points out it isn't already there.

[clamus]

simulation.sources has 2 ways to simulate data. One is simulate_sparse_stc, where a either a single dipole is placed in each label, or several dipoles are randomly places, and the same "time series" is set for all dipoles. The other is simulate_stc, where one gives the time course in each label and a spatial weighting kernel to have like a distributed source in each label. Either way, there is a kind of restriction of the spatial and temporal traces that can be simulated. What I need is no restrictions at all. I want to select areas or cortex and place on those "distributed" areas arbitrary time series, with arbitrary spatial and temporal correlations.

[larsoner]

Ahh you're right, I thought we had more general support. It wouldn't be too bad to add to simulate_stc or make a new function that allowed you to pass a SourceEstimate and fwd to combine to simulate data. What that take care of your use case?

[clamus]

Yes, that would definitely do it! I can take care of it is this is OK. Also, I finish fixing the issue with restrict_forward_to_label

[dgwakeman]

For what it's worth is it even worth generating a function for that. Just multiply the two matrices and bam you have your sensor data, slam it into a template file and you are done. it should be like 3 lines of code.

[larsoner]

I agree the function should be pretty simple -- do you volunteer to make the PR, then? :)

[clamus]

Yup, I volunteer :) And agree this is pretty simple.