Hi all,
I am trying to solve two related problems during source localization. First, I’m trying to extract all vertices in a given parcel of a specific parcellation, convert into MNI coordinates, and then, I’m trying to find the centroid of the parcel. Based on code from another response on this forum, I have managed to extract all vertices in a given parcel, albeit I have some concerns whether I am doing that correctly. Please see the attached code below - could someone confirm whether version 1, 2, 3 or 4 is the appropriate approach?
Finally, I know that the way I am calculating the centroid of the parcel is incorrect - I’m simply taking the mean of all coordinates, even though the surface is curved, not flat. However, I am unsure how to find the centroid in a situation such as this. Could someone help? Thank you so much in advance.
# Set up source space
# sub here is the subject ID, and subjects_dir is the MRI folder
src = mne.setup_source_space(
subject=sub, subjects_dir=subjects_dir,
spacing='oct6', add_dist=False)
labels = mne.read_labels_from_annot(sub,parc ='parc2018yeo17_200',subjects_dir = subjects_dir)
labelname = '17Networks_' + hemifield.upper() + '_' + roi + '-' + hemifield
offset = src[0]['vertno'].shape[0] # this value will be used below, and is based on
# code from this link:
#https://mne.discourse.group/t/correspondence-of-sources-with-parcellation/1942/5
for label in labels:
if label.name == labelname:
#print(label)
this_hemi = 0 if (label.hemi == 'lh') else 1
# identify vertices that are present in the src of the participant
# AND belong to the vertices in a given parcellation based on the
# label information:
idx = np.intersect1d(label.vertices, src[this_hemi]['vertno'])
# not sure about why searchsorted is necessary:
parc_hemi_idx = np.searchsorted(src[this_hemi]['vertno'], idx)
# add offset if in right hemisphere:
parc_idx = parc_hemi_idx + offset * this_hemi
# different ways to get MNI coordinates:
# 1) use the idx without the hemi specific index and simply specify the hemi
# through the second argument in vertex_to_mni (this should only matter for rh)
mnicoordA = mne.vertex_to_mni(parc_hemi_idx,this_hemi,sub,subjects_dir=subjects_dir)
# 2) use the vertex info from the label, without reference to the src
# I guesss this is problematic, because the label might have vertices that the src does not?
mnicoordB = mne.vertex_to_mni(label.vertices,this_hemi,sub,subjects_dir=subjects_dir)
# 3) use the idx of the vertices that overlap between the label and the src directly
mnicoordC = mne.vertex_to_mni(idx,this_hemi,sub,subjects_dir=subjects_dir)
# 4) index the src with the searchsorted idx, and look up these vertices
mnicoordD = mne.vertex_to_mni(src[this_hemi]['vertno'][parc_hemi_idx],this_hemi,sub,subjects_dir=subjects_dir)
# 3 and 4 seem to produce identical results so I assume they do the same thing
# centroids coming from the diff methods:
centroidA = np.mean(mnicoordA,axis=0)#
centroidB = np.mean(mnicoordB,axis=0)
centroidC = np.mean(mnicoordC,axis=0)
centroidD = np.mean(mnicoordD,axis=0)
# A seems way off, C and D are identical, B is comparable to C and D
break