Extracting MNI coordinates of vertices in a specific parcel

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
  • MNE version:1.1.1
  • operating system: Windows 10

I have found this built-in function, could someone confirm that these two lines are indeed doing what I hope they are doing, i.e., identifying (roughly) the center of a parcel, and converting that point to MNI coordinates?

mne_vert = label.center_of_mass(subject=sub,restrict_vertices=src,subjects_dir=subjects_dir)
centroidX = mne.vertex_to_mni(mne_vert,this_hemi,sub,subjects_dir=subjects_dir)

those lines look correct to me. resulting centroidX should be a 1x3 array of XYZ in millimeters.

Thanks, Dan. Indeed, centroidX is an array of coordinates. Is my interpretation then correct that these lines will get the center of mass of all vertices within a label that are also present in the src, and then transform this vertex into regular old MNI coordinates (in mm)? Will mne_vert in my example also always be part of the vertices in src?

mne_vert should be guaranteed to be one of the vertices in src, yes. But your description of how it’s found is I think a little inaccurate:

  • I think you said “it will find all points in the label that are also in src and then compute the centroid of those points”
  • I think the correct explanation is “it will calculate the centroid based only on the vertices in the label and then pick the point in your src that is closest to that centroid”

but I’d have to double-check the code to be absolutely sure.

Thanks, that makes perfect sense too.