split_parc.py

from pathlib import Path
import numpy as np
import mne
import nibabel as nib
import matplotlib.pyplot as plt

def split_parc(xyz, label, n, axis='y'):

    axes = dict(x=0, y=1, z=2)

    m = xyz[:, axes[axis]].min()
    M = xyz[:, axes[axis]].max()
    bounds = (M-m) * np.linspace(0, 1., 1+n) + m

    groups = np.digitize(xyz[:, axes[axis]], 
                         bounds)
    
    labels = list()
    for group_id, _ in enumerate(bounds):
        label_ = label.copy()
        label_.name = area + f'_{group_id}'
        label_.vertices = label.vertices[groups==group_id]
        labels.append(label_)
    return labels

if __name__ == '__main__':

    SUBJECTS_DIR = Path('/Users/jeanremi/git/plot_brain/plot_brain/data/')
    hemi = 'lh'
    surf = 'pial'
    parc = 'PALS_B12_Brodmann'

    surf = SUBJECTS_DIR / 'fsaverage' / 'surf' / f'{hemi}.{surf}'

    # read geometry
    xyz, tris = nib.freesurfer.read_geometry(surf)

    # read parcellation
    labels = mne.read_labels_from_annot('fsaverage', 
                                        parc=parc, 
                                        subjects_dir=SUBJECTS_DIR, 
                                        verbose=False)

    # select area
    area = f'Brodmann.22-{hemi}'
    label = [l for l in labels if l.name == area][0]

    split_label = split_parc(xyz[label.vertices], label, 3)
    
    # plot
    plt.scatter(xyz[::100, 1],
            xyz[::100, 2], c='gray')
    for l in split_label:
        plt.scatter(xyz[l.vertices, 1],
                    xyz[l.vertices, 2])