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.
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
@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
@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.
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.
These typically do not do plan averages in the parcel. One key reason is that sources have a directions so by averaging the raw time courses obtained with fixed orientation or pick normal you have cancellation effects. You will have the same problem with your approach as you have columns of the forward solution and these can be basically of opposite polarity for dipoles eg located on opposite walls of a gyrus.