Compute adjacency between MNE Label objects.

label_adjacency.py

import matplotlib.pyplot as plt
import mne
import numpy as np
from scipy import sparse

# We will operate on the fsaverage source space. Fetch it first.
subjects_dir = mne.datasets.sample.data_path() / "subjects"
mne.datasets.fetch_fsaverage(subjects_dir)

# Load the source space.
src = mne.read_source_spaces(
    subjects_dir / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
)

# Fetch and load the HCP labels
mne.datasets.fetch_hcp_mmp_parcellation(subjects_dir=subjects_dir, verbose=True)
labels = mne.read_labels_from_annot(
    "fsaverage", "HCPMMP1", "lh", subjects_dir=subjects_dir
)


def label_adjacency(labels, src):
    """Compute adjacency between labels.

    Two labels are considered adjacent if one of their vertices are adjacent.

    Parameters
    ----------
    labels : list of mne.Label
        The labels between which to compute adjacency.
    src : mne.SourceSpaces
        The source space on which the labels are defined.

    Returns
    -------
    label_adjacency : scipy.sparse.coo_matrix
        A sparse adjacency matrix containing a 1 for labels that are adjacenct and 0
        otherwise.
    """
    src_adjacency = mne.spatial_src_adjacency(src).tocsr()
    label_verts = list()
    for label in labels:
        src_hemi = src[0] if label.hemi == "lh" else src[1]
        label_verts.append(label.get_vertices_used(src_hemi["vertno"]))

    adjacent_label_inds = list()  # list of pairs of label indices
    for ind1, label1 in enumerate(labels):
        for ind2, label2 in enumerate(labels):
            # If the labels are on different hemispheres, they are not adjacent.
            if label1.hemi != label2.hemi:
                continue

            # Get adjacent vertices if any.
            adj_verts = src_adjacency[label_verts[ind1]][:, label_verts[ind2]]
            if adj_verts.count_nonzero() > 0:
                adjacent_label_inds.append((ind1, ind2))
    return sparse.coo_matrix(
        (np.ones(len(adjacent_label_inds)), tuple(zip(*adjacent_label_inds))),
        shape=(len(labels), len(labels)),
    )


# Use the function above to compute adjacency between labels.
adjacency = label_adjacency(labels, src)

# Visualize the adjacency.
plt.imshow(adjacency.todense())