LeadField matrix

Hi everyone, I am trying to obtain a leadfield matrix (channels x brain regions) which in my case would be 62x200 (I am using a brain atlas with 200 parcels). Any idea on how to move from the “canonical” leadfield matrix (channels x sources) to the one I’ve just described. Since mne has a function to extract the parcels time series, there should be an easy way to extract the leadfield matrix without scripting it from scratch.

Thanks in advance

what would be for you the forward field of a parcel? putting the full parcel with a constant activation of 1?

Alex

Hi @agramfort thanks for reply to this. It would be the same as the “canonical” leadfield array but instead having the vertices on the second dim, it will have the average of the values for the vertices that falls in the parcels. Anyway, I’ve just solved the issue myself…just some extra-work with mne.decimate

Thanks

ok that’s equivalent to what I suggest above. It’s as if you were to put a constant activity over each parcel.

you would just want the leadfield as a numpy array?

Alex

@agramfort please find the piece of code I wrote below for a sanity check and for other people who might encounter the same issue. Could you please take a look at it?

stcs = mne.read_source_estimate(op.join(eeg_path + 'sub_0014-stc'))

lh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/lh.Schaefer2018_200Parcels_7Networks_order.annot'))[0]

rh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/rh.Schaefer2018_200Parcels_7Networks_order.annot'))[0]

regmap = np.concatenate([lh_regmap, rh_regmap+np.max(lh_regmap)])

vtx2use = np.concatenate([stcs.vertices[0], stcs.vertices[1]+np.max(stcs.vertices[0])])
regions_idx = regmap[vtx2use]


fwd = mne.read_forward_solution(op.join(eeg_path + 'sub_0014-fwd'))
leadfield = fwd['sol']['data']
print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape)



fwd_fixed = mne.convert_forward_solution(fwd, surf_ori=True, force_fixed=True,
                                         use_cps=True)
leadfield = fwd_fixed['sol']['data']
print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape)



new_leadfield = []
for xx in range(1,np.unique(regions_idx).shape[0]):
    new_leadfield.append(np.mean(np.squeeze(leadfield[:,np.where(regions_idx==xx)]),axis=1))
    
new_leadfield = np.array(new_leadfield).T

print("Leadfield size : %d sensors x %d regions" % new_leadfield.shape)```


Output is: 

Leadfield size : 62 sensors x 23121 dipoles
Leadfield size : 62 sensors x 7707 dipoles
Leadfield size : 62 sensors x 200 regions


Could you please confirm that is correct? 

Thanks

it’s hard to tell without running the code. Can you share snippet using the mne-sample-data so it works also on my machine?

thanks
Alex

@agramfort I apologize but I don’t have the mne-sample-data here. Anyway it should be easy to run. Please set up:
stcs = the source estimate
lh_regmap = path to the annotation file of the .annot file (in your case it would be the lh.aparc.annot)
rh_regmap = same as lh_regmap but for the right hemisphere.

Here is the code.

stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method, pick_ori=None, return_generator=True)
lh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/lh.aparc.annot'))[0]

rh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/rh.aparc.annot'))[0]

Hope this helps

hi Davide,

yes it’s easy for me to do but it’s still minutes for me. You can easily do this yourself too.
Let’s balance the efforts here so we can collectively progress.

Alex

1 Like

You can download the data by running mne.datasets.sample.data_path()