diff --git a/doc/api/statistics.rst b/doc/api/statistics.rst
index c6f174152b3..90eca48d7f1 100644
--- a/doc/api/statistics.rst
+++ b/doc/api/statistics.rst
@@ -68,6 +68,7 @@ Compute ``adjacency`` matrices for cluster-level statistics:
 
    channels.find_ch_adjacency
    channels.read_ch_adjacency
+   label_adjacency
    spatial_dist_adjacency
    spatial_src_adjacency
    spatial_tris_adjacency
@@ -75,3 +76,4 @@ Compute ``adjacency`` matrices for cluster-level statistics:
    spatio_temporal_src_adjacency
    spatio_temporal_tris_adjacency
    spatio_temporal_dist_adjacency
+   volume_label_adjacency
diff --git a/doc/changes/dev/14226.newfeature.rst b/doc/changes/dev/14226.newfeature.rst
new file mode 100644
index 00000000000..dfe104afcdb
--- /dev/null
+++ b/doc/changes/dev/14226.newfeature.rst
@@ -0,0 +1,3 @@
+Add :func:`mne.label_adjacency` and :func:`mne.volume_label_adjacency` for determining
+whether labels are adjacent (touching) or not, by :newcontrib:`Raphaël Bordas` and `Marijn van
+Vliet`_
diff --git a/doc/changes/names.inc b/doc/changes/names.inc
index 6fd9a97073b..6856ca7e2da 100644
--- a/doc/changes/names.inc
+++ b/doc/changes/names.inc
@@ -364,6 +364,7 @@
 .. _Ram Pari: https://github.com/ramkpari
 .. _Ramiro Gatti: https://github.com/ragatti
 .. _ramonapariciog: https://github.com/ramonapariciog
+.. _Raphaël Bordas: https://github.com/raphbrd
 .. _Rasmus Aagaard: https://github.com/rasgaard
 .. _Rasmus Zetter: https://github.com/rzetter
 .. _Reza Nasri: https://github.com/rznas
diff --git a/mne/__init__.pyi b/mne/__init__.pyi
index a7ea80344d7..1b5821536a5 100644
--- a/mne/__init__.pyi
+++ b/mne/__init__.pyi
@@ -90,6 +90,7 @@ __all__ = [
     "head_to_mri",
     "inverse_sparse",
     "io",
+    "label_adjacency",
     "label_sign_flip",
     "labels_to_stc",
     "make_ad_hoc_cov",
@@ -188,6 +189,7 @@ __all__ = [
     "verbose",
     "vertex_to_mni",
     "viz",
+    "volume_label_adjacency",
     "what",
     "whiten_evoked",
     "write_bem_solution",
@@ -351,6 +353,7 @@ from .label import (
     BiHemiLabel,
     Label,
     grow_labels,
+    label_adjacency,
     label_sign_flip,
     labels_to_stc,
     morph_labels,
@@ -359,6 +362,7 @@ from .label import (
     read_labels_from_annot,
     split_label,
     stc_to_label,
+    volume_label_adjacency,
     write_label,
     write_labels_to_annot,
 )
diff --git a/mne/label.py b/mne/label.py
index 6907fac05c2..3eb921ba6be 100644
--- a/mne/label.py
+++ b/mne/label.py
@@ -8,10 +8,12 @@
 import re
 from collections import defaultdict
 from colorsys import hsv_to_rgb, rgb_to_hsv
+from pathlib import Path
 
 import numpy as np
 from scipy import linalg
 
+from ._freesurfer import get_volume_labels_from_aseg
 from .fixes import _safe_svd
 from .morph_map import read_morph_map
 from .parallel import parallel_func
@@ -19,6 +21,7 @@
     SourceEstimate,
     VolSourceEstimate,
     _center_of_mass,
+    _volume_labels,
     extract_label_time_course,
     spatial_src_adjacency,
 )
@@ -3043,3 +3046,130 @@ def select_sources(
         )
 
     return new_label
+
+
+def label_adjacency(labels, src):
+    """Compute adjacency between labels.
+
+    Two labels are considered adjacent if one of their vertices are adjacent in the
+    source space.
+
+    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 adjacent and 0
+        otherwise.
+
+    See Also
+    --------
+    volume_label_adjacency
+
+    Notes
+    -----
+    .. versionadded:: 1.13
+    """
+    from scipy.sparse import coo_matrix
+
+    src_adjacency = spatial_src_adjacency(src).tocsr()
+    label_src_ind = list()
+    for label in labels:
+        src_hemi = src[0] if label.hemi == "lh" else src[1]
+        label_verts = label.get_vertices_used(src_hemi["vertno"])
+        src_ind = np.searchsorted(src_hemi["vertno"], label_verts)
+        if label.hemi == "rh":
+            src_ind += src[0]["nuse"]
+        label_src_ind.append(src_ind)
+
+    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_src_ind[ind1], :][:, label_src_ind[ind2]]
+            if adj_verts.data.any():
+                adjacent_label_inds.append((ind1, ind2))
+    return coo_matrix(
+        (np.ones(len(adjacent_label_inds)), tuple(zip(*adjacent_label_inds))),
+        shape=(len(labels), len(labels)),
+    )
+
+
+@fill_doc
+def volume_label_adjacency(src, subject, subjects_dir, *, aseg="auto", labels=None):
+    """Compute adjacency between volume labels.
+
+    Two labels are considered adjacent if one of their voxels are adjacent in the
+    (volumetric) source space.
+
+    Parameters
+    ----------
+    src : mne.SourceSpaces
+        The volumetric source space on which the labels are defined.
+    %(subject)s
+    %(subjects_dir)s
+    %(aseg)s
+    %(labels_aseg)s
+
+    Returns
+    -------
+    label_adjacency : scipy.sparse.coo_matrix
+        A sparse adjacency matrix containing a 1 for labels that are adjacent and 0
+        otherwise.
+    labels : list of str
+        The names of the labels which contain at least one source point.
+
+    See Also
+    --------
+    label_adjacency
+
+    Notes
+    -----
+    .. versionadded:: 1.13
+    """
+    from scipy import sparse
+
+    subjects_dir = Path(get_subjects_dir(subjects_dir, raise_error=True))
+    if aseg == "auto":  # use aparc+aseg if auto
+        aseg = _check_fname(
+            subjects_dir / subject / "mri" / "aparc+aseg.mgz",
+            overwrite="read",
+            must_exist=False,
+        )
+        if not aseg:  # if doesn't exist use wmparc
+            aseg = subjects_dir / subject / "mri" / "wmparc.mgz"
+    else:
+        aseg = subjects_dir / subject / "mri" / f"{aseg}.mgz"
+
+    if labels is None:
+        labels = get_volume_labels_from_aseg(aseg)
+
+    vol_labels = _volume_labels(src, (aseg, labels), mri_resolution=False)
+    src_adjacency = spatial_src_adjacency(src).tocsr()
+
+    label_verts = list()
+    for label in vol_labels:
+        label_verts.append(np.searchsorted(src[0]["vertno"], label.vertices))
+
+    adjacent_label_inds = list()  # list of pairs of label indices
+    for ind1, verts1 in enumerate(label_verts):
+        for ind2, verts2 in enumerate(label_verts):
+            # Get adjacent vertices if any.
+            adj_verts = src_adjacency[verts1, :][:, verts2]
+            if adj_verts.data.any():
+                adjacent_label_inds.append((ind1, ind2))
+
+    adj = sparse.coo_matrix(
+        (np.ones(len(adjacent_label_inds)), tuple(zip(*adjacent_label_inds))),
+        shape=(len(labels), len(labels)),
+    )
+    return adj, labels
diff --git a/mne/tests/test_label.py b/mne/tests/test_label.py
index eca3cc5ee07..b485b4a476f 100644
--- a/mne/tests/test_label.py
+++ b/mne/tests/test_label.py
@@ -20,6 +20,7 @@
 )
 from scipy import sparse
 
+import mne
 from mne import (
     grow_labels,
     labels_to_stc,
@@ -57,6 +58,9 @@
 real_label_fname = data_path / "MEG" / "sample" / "labels" / "Aud-lh.label"
 v1_label_fname = subjects_dir / "sample" / "label" / "lh.V1.label"
 
+fname_vsrc = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-vol-7-fwd.fif"
+fname_src_fs = data_path / "subjects" / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
+
 fwd_fname = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-6-fwd.fif"
 src_bad_fname = data_path / "subjects" / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
 label_dir = subjects_dir / "sample" / "label" / "aparc"
@@ -64,6 +68,7 @@
 test_path = Path(__file__).parents[1] / "io" / "tests" / "data"
 label_fname = test_path / "test-lh.label"
 
+
 # This code was used to generate the "fake" test labels:
 # for hemi in ['lh', 'rh']:
 #    label = Label(np.unique((np.random.rand(100) * 10242).astype(int)),
@@ -1263,3 +1268,119 @@ def test_label_geometry(fname, area):
     )
     assert_array_less(inside_euc, inside_dist)
     assert_array_less(0.25 * inside_dist, inside_euc)
+
+
+@testing.requires_testing_data
+def test_volume_label_adjacency():
+    """Test label adjacency."""
+    pytest.importorskip("nibabel")
+    pytest.importorskip("sklearn")
+    src = read_source_spaces(fname_vsrc)
+
+    # aseg=auto uses the aparc+aseg atlas, which does not exist in the testing datasets
+    adj, labels = mne.volume_label_adjacency(
+        src, subject="sample", subjects_dir=subjects_dir, aseg="aseg"
+    )
+    n_neighbors = adj.sum(axis=1)
+
+    assert_equal(len(labels), 46)  # default number of labels in aseg.mgz
+    assert_equal(adj.shape, (len(labels), len(labels)))
+
+    assert_equal(n_neighbors.min(), 0)
+    assert_equal(np.sum(n_neighbors == 0), 4)
+
+    # example: 'Left-Thalamus-Proper'
+    label_idx = 7
+    connected_labels_idx = adj.toarray()[label_idx, :]
+    connected_labels = np.array(labels)[np.where(connected_labels_idx == 1)[0]]
+
+    assert_equal(
+        np.sort(connected_labels).tolist(),
+        [
+            "3rd-Ventricle",
+            "Brain-Stem",
+            "CSF",
+            "Left-Accumbens-area",
+            "Left-Cerebral-Cortex",
+            "Left-Cerebral-White-Matter",
+            "Left-Hippocampus",
+            "Left-Lateral-Ventricle",
+            "Left-Thalamus-Proper",
+            "Left-VentralDC",
+            "Unknown",
+        ],
+    )
+
+    input_labels = [
+        "Left-Thalamus-Proper",
+        "Left-Hippocampus",
+        "Right-Hippocampus",
+    ]
+    adj, labels = mne.volume_label_adjacency(
+        src,
+        subject="sample",
+        subjects_dir=subjects_dir,
+        aseg="aseg",
+        labels=input_labels,
+    )
+
+    assert_equal(
+        adj.toarray(),
+        np.array(
+            [
+                [1, 1, 0],
+                [1, 1, 0],
+                [0, 0, 1],
+            ]
+        ),
+    )
+
+    assert_equal(labels, input_labels)
+
+    with pytest.raises(FileNotFoundError):
+        mne.volume_label_adjacency(
+            src,
+            subject="sample",
+            subjects_dir=subjects_dir,
+            aseg="my-aseg",
+            labels=input_labels,
+        )
+
+
+@testing.requires_testing_data
+def test_label_adjacency():
+    """Test label adjacency."""
+    pytest.importorskip("nibabel")
+    src = read_source_spaces(fname_src_fs)
+    mne.add_source_space_distances(src, dist_limit=0.01, n_jobs=-1)
+
+    labels = mne.read_labels_from_annot(
+        subject="fsaverage",
+        subjects_dir=subjects_dir,
+    )
+    adj = mne.label_adjacency(labels, src)
+
+    n_neighbors = adj.sum(axis=1)
+
+    assert_equal(len(labels), 69)  # default number of labels in aseg.mgz
+    assert_equal(adj.shape, (len(labels), len(labels)))
+
+    assert_equal(n_neighbors.min(), 0)
+    assert_equal(np.sum(n_neighbors == 0), 1)
+
+    input_labels = [
+        "cuneus-lh",
+        "cuneus-rh",
+        "precuneus-lh",
+    ]
+    adj = mne.label_adjacency([lab for lab in labels if lab.name in input_labels], src)
+    assert_equal(
+        adj.toarray(),
+        np.array(
+            [
+                [1, 0, 1],
+                [0, 1, 0],
+                [1, 0, 1],
+            ]
+        ),
+    )
diff --git a/mne/utils/docs.py b/mne/utils/docs.py
index 8ae9ff682ba..7143b6a1d28 100644
--- a/mne/utils/docs.py
+++ b/mne/utils/docs.py
@@ -2497,7 +2497,13 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
 label_tc : array | list (or generator) of array, shape (n_labels[, n_orient], n_times)
     Extracted time course for each label and source estimate.
 """
-
+docdict["labels_aseg"] = """
+labels : list of str | None
+    Labeled regions of interest to plot. See :func:`mne.get_montage_volume_labels`
+    for one way to determine regions of interest. Regions can also be chosen from
+    the :term:`FreeSurfer LUT`. If ``None``, all labels that are defined in the
+    segmentation file are used.
+"""
 docdict["labels_eltc"] = """
 labels : Label | BiHemiLabel | list | tuple | str
     If using a surface or mixed source space, this should be the
diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py
index bbbd6bc550c..b296fc0a446 100644
--- a/mne/viz/_brain/_brain.py
+++ b/mne/viz/_brain/_brain.py
@@ -2935,11 +2935,7 @@ def add_volume_labels(
         Parameters
         ----------
         %(aseg)s
-        labels : list
-            Labeled regions of interest to plot. See
-            :func:`mne.get_montage_volume_labels`
-            for one way to determine regions of interest. Regions can also be
-            chosen from the :term:`FreeSurfer LUT`.
+        %(labels_aseg)s
         colors : list | matplotlib-style color | None
             A list of anything matplotlib accepts: string, RGB, hex, etc.
             (default :term:`FreeSurfer LUT` colors).