Deep source modeling (phase 1)

doc/api/forward.rst

@@ -39,6 +39,7 @@ Forward Modeling
    read_surface
    sensitivity_map
    setup_source_space
+   setup_subcortical_source_space
    setup_volume_source_space
    surface.complete_surface_info
    surface.read_curvature

mne/__init__.pyi

@@ -164,6 +164,7 @@ __all__ = [
     "set_log_level",
     "set_memmap_min_size",
     "setup_source_space",
+    "setup_subcortical_source_space",
     "setup_volume_source_space",
     "simulation",
     "source_space",
@@ -404,6 +405,7 @@ from .source_space._source_space import (
     morph_source_spaces,
     read_source_spaces,
     setup_source_space,
+    setup_subcortical_source_space,
     setup_volume_source_space,
     write_source_spaces,
 )

mne/_fiff/constants.py

@@ -401,6 +401,7 @@
 FIFF.FIFFV_MNE_SURF_LEFT_HEMI = 101
 FIFF.FIFFV_MNE_SURF_RIGHT_HEMI = 102
 FIFF.FIFFV_MNE_SURF_MEG_HELMET = 201  # Use this irrespective of the system
+FIFF.FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE = 1000  # + aseg value, e.g. hippocampus
 #
 #   These relate to the Isotrak data (enum(point))
 #

mne/_fiff/tests/test_constants.py

@@ -57,7 +57,10 @@
     "viewkeys",
     "viewvalues",  # Py2
 )
-_tag_ignore_names = ()  # for fiff-constants pending updates
+_tag_ignore_names = (
+    # pending addition to fiff-constants, see mne-tools/mne-python#14130
+    "FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE",
+)
 _ignore_incomplete_enums = (  # XXX eventually we could complete these
     "bem_surf_id",
     "cardinal_point_cardiac",

mne/source_estimate.py

@@ -415,7 +415,7 @@ def _get_src_type(src, vertices, warn_text=None):
             src_type = "mixed"
     else:
         src_type = src.kind
-    assert src_type in ("surface", "volume", "mixed", "discrete")
+    assert src_type in ("surface", "volume", "mixed", "discrete", "subcortical_surf")
     return src_type
 
 
@@ -443,7 +443,7 @@ def guess_src_type():
     # infer Klass from src_type
     if src_type == "surface":
         Klass = VectorSourceEstimate if vector else SourceEstimate
-    elif src_type in ("volume", "discrete"):
+    elif src_type in ("volume", "discrete", "subcortical_surf"):
         Klass = VolVectorSourceEstimate if vector else VolSourceEstimate
     elif src_type == "mixed":
         Klass = MixedVectorSourceEstimate if vector else MixedSourceEstimate

mne/source_space/__init__.pyi

@@ -6,6 +6,7 @@ __all__ = [
     "get_decimated_surfaces",
     "read_source_spaces",
     "setup_source_space",
+    "setup_subcortical_source_space",
     "setup_volume_source_space",
     "write_source_spaces",
 ]
@@ -17,6 +18,7 @@ from ._source_space import (
     get_decimated_surfaces,
     read_source_spaces,
     setup_source_space,
+    setup_subcortical_source_space,
     setup_volume_source_space,
     write_source_spaces,
 )

mne/source_space/_source_space.py

@@ -35,6 +35,7 @@
 )
 from .._freesurfer import (
     _check_mri,
+    _get_aseg,
     _get_atlas_values,
     _get_mri_info_data,
     get_volume_labels_from_aseg,
@@ -50,6 +51,8 @@
     _create_surf_spacing,
     _get_ico_surface,
     _get_surf_neighbors,
+    _keep_largest_component,
+    _marching_cubes,
     _normalize_vectors,
     _tessellate_sphere_surf,
     _triangle_neighbors,
@@ -299,6 +302,7 @@ def __init__(self, source_spaces, info=None):
     @property
     def kind(self):
         types = list()
+        ids = list()
         for si, s in enumerate(self):
             _validate_type(s, dict, f"source_spaces[{si}]")
             types.append(s.get("type", None))
@@ -307,19 +311,40 @@ def kind(self):
                 types[-1],
                 ("surf", "discrete", "vol"),
             )
-        if all(k == "surf" for k in types[:2]):
+            ids.append(s.get("id", FIFF.FIFFV_MNE_SURF_UNKNOWN))
+        n = len(types)
+        is_subcortical = [
+            t == "surf" and i >= FIFF.FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE
+            for t, i in zip(types, ids)
+        ]
+        leading_surf_pair = n >= 2 and types[0] == "surf" and types[1] == "surf"
+        leading_subcortical_pair = (
+            leading_surf_pair and is_subcortical[0] and is_subcortical[1]
+        )
+        if leading_surf_pair and not leading_subcortical_pair:
             surf_check = 2
-            if len(types) == 2:
-                kind = "surface"
-            else:
-                kind = "mixed"
+            kind = "surface" if n == 2 else "mixed"
+        elif n == 1 and types[0] == "surf" and not is_subcortical[0]:
+            surf_check = 1
+            kind = "mixed"
+        elif n == 0:
+            surf_check = 0
+            kind = "mixed"
+        elif all(is_subcortical):
+            surf_check = 0
+            kind = "subcortical_surf"
+        elif any(is_subcortical):
+            surf_check = 0
+            kind = "mixed"
+        elif all(k == "discrete" for k in types):
+            surf_check = 0
+            kind = "discrete"
         else:
             surf_check = 0
-            if all(k == "discrete" for k in types):
-                kind = "discrete"
-            else:
-                kind = "volume"
-        if any(k == "surf" for k in types[surf_check:]):
+            kind = "volume"
+        if any(
+            types[i] == "surf" and not is_subcortical[i] for i in range(surf_check, n)
+        ):
             raise RuntimeError(f"Invalid source space with kinds {types}")
         return kind
 
@@ -1036,6 +1061,7 @@ def _read_one_source_space(fid, this):
                 offset += n
             res["neighbor_vert"] = neighbors
 
+    if res["type"] in ("vol", "surf"):
         tag = find_tag(fid, this, FIFF.FIFF_COMMENT)
         if tag is not None:
             res["seg_name"] = tag.data
@@ -1432,7 +1458,7 @@ def _write_one_source_space(fid, this, verbose=None):
         )
 
     #   Segmentation data
-    if this["type"] == "vol" and ("seg_name" in this):
+    if this["type"] in ("vol", "surf") and ("seg_name" in this):
         # Save the name of the segment
         write_string(fid, FIFF.FIFF_COMMENT, this["seg_name"])
 
@@ -1982,6 +2008,175 @@ def _complete_vol_src(sp, subject=None):
     return sp
 
 
+def _surf_from_mesh(rr, tris, subject):
+    """Build a source-space-ready surf dict from vertices/triangles (in m)."""
+    surf = dict(rr=np.asarray(rr, float), tris=np.asarray(tris, np.int64))
+    complete_surface_info(surf, do_neighbor_vert=False, copy=False)
+    surf["inuse"] = np.ones(surf["np"], int)
+    sizes = _normalize_vectors(surf["nn"])
+    surf["inuse"][sizes <= 0] = False
+    surf["nuse"] = int(surf["inuse"].sum())
+    surf["vertno"] = np.where(surf["inuse"])[0]
+    surf["use_tris"] = None
+    surf["nuse_tri"] = 0
+    surf["subject_his_id"] = subject
+    for key in ("tri_area", "tri_cent", "tri_nn", "neighbor_tri"):
+        del surf[key]
+    surf.update(
+        dist=None,
+        dist_limit=None,
+        nearest=None,
+        nearest_dist=None,
+        pinfo=None,
+        patch_inds=None,
+        type="surf",
+        coord_frame=FIFF.FIFFV_COORD_MRI,
+    )
+    return surf
+
+
+@verbose
+def setup_subcortical_source_space(
+    subject,
+    label=None,
+    surface=None,
+    aseg="auto",
+    subjects_dir=None,
+    keep_largest_component=True,
+    smooth=0,
+    fill_hole_size=None,
+    add_dist=False,
+    *,
+    verbose=None,
+):
+    """Set up a subcortical or cerebellar surface source space.
+
+    This builds a :class:`~mne.SourceSpaces` from a triangulated mesh of a subcortical
+    or cerebellar structure, either tessellated directly from an
+    anatomical segmentation (``label``) or supplied as an
+    externally-produced mesh (``surface``, e.g. one fitted by another
+    package such as CMB). Exactly one of ``label`` or ``surface`` must be
+    provided.
+
+    Parameters
+    ----------
+    subject : str
+        Subject to process.
+    label : str | list | dict | None
+        Region(s) of interest to tessellate from the anatomical
+        segmentation given by ``aseg``. One source space is created per
+        entry (a single str is turned into a one-element list). If dict,
+        maps region names to atlas id numbers, allowing the use of other
+        atlases. Mutually exclusive with ``surface``.
+    surface : path-like | dict | None
+        A FreeSurfer-compatible surface file (e.g. a ``.surf`` file), or a
+        dict with ``'rr'`` and ``'tris'`` entries in FreeSurfer surface RAS
+        coordinates (mm), such as those returned by :func:`mne.read_surface`
+        or produced by an external mesh-fitting tool. Creates a single
+        source space. Mutually exclusive with ``label``.
+    %(aseg)s
+        Only used when ``label`` is provided.
+    %(subjects_dir)s
+    keep_largest_component : bool
+        If True (default), keep only the largest connected component of
+        each tessellated mesh, discarding disconnected islands (the
+        marching-cubes equivalent of FreeSurfer's
+        ``mris_extract_main_component``).
+    %(smooth)s
+        Only used when ``label`` is provided.
+    fill_hole_size : int | None
+        The size of holes to remove in the mesh in voxels. Default is None,
+        no holes are removed. This dilates the boundaries of the surface by
+        ``fill_hole_size`` voxels, so use the minimal size needed. Only used
+        when ``label`` is provided.
+    add_dist : bool
+        If True, compute inter-source distances along the mesh (see
+        :func:`mne.add_source_space_distances`). Default False, as this can
+        be slow and is not needed for a forward solution.
+    %(verbose)s
+
+    Returns
+    -------
+    src : instance of SourceSpaces
+        The subcortical/cerebellar surface source space(s), one per
+        ``label`` entry, or a single one if ``surface`` was used.
+
+    See Also
+    --------
+    setup_volume_source_space
+    setup_source_space
+
+    Notes
+    -----
+    This is a first, deliberately narrow proof of concept: it has been
+    validated interactively on the ``sample`` subject. Known gaps, to be
+    addressed in follow-up work: morphing (:class:`~mne.SourceMorph`) does
+    not yet support these source spaces, and
+    :func:`mne.extract_label_time_course` does not yet know how to select
+    vertices within a subcortical-surface label.
+
+    .. versionadded:: 1.12
+    """
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    _validate_type(label, (str, list, tuple, dict, None), "label")
+    _validate_type(surface, ("path-like", dict, None), "surface")
+    if (label is None) == (surface is None):
+        raise ValueError(
+            "Exactly one of `label` or `surface` must be provided, got "
+            f"label={label!r}, surface={surface!r}"
+        )
+
+    srcs = list()
+    if label is not None:
+        aseg_img, aseg_data = _get_aseg(aseg, subject, subjects_dir)
+        mri = aseg_img.get_filename()
+        volume_label = _check_volume_labels(label, mri, name="label")
+        vox_mri_t = np.array(aseg_img.header.get_vox2ras_tkr(), float)
+        vox_mri_t[:3] *= 1e-3  # mm -> m
+        meshes = _marching_cubes(
+            aseg_data,
+            list(volume_label.values()),
+            smooth=smooth,
+            fill_hole_size=fill_hole_size,
+        )
+        for (seg_name, seg_id), (rr, tris) in zip(volume_label.items(), meshes):
+            if len(rr) == 0:
+                warn(
+                    f"Value {seg_id} not found for label {seg_name!r} in "
+                    f"anatomical segmentation file {mri}, skipping"
+                )
+                continue
+            if keep_largest_component:
+                rr, tris = _keep_largest_component(rr, tris)
+            rr = apply_trans(vox_mri_t, rr)
+            s = _surf_from_mesh(rr, tris, subject)
+            s["seg_name"] = seg_name
+            s["id"] = FIFF.FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE + seg_id
+            srcs.append(s)
+        if len(srcs) == 0:
+            raise ValueError(f"None of the requested labels were found in {mri}")
+    else:
+        if isinstance(surface, dict):
+            rr, tris = surface["rr"], surface["tris"]
+        else:
+            surface = str(
+                _check_fname(surface, overwrite="read", must_exist=True, name="surface")
+            )
+            rr, tris = read_surface(surface)[:2]
+        rr = np.array(rr, float) / 1000.0  # mm -> m
+        tris = np.array(tris, np.int64)
+        if keep_largest_component:
+            rr, tris = _keep_largest_component(rr, tris)
+        s = _surf_from_mesh(rr, tris, subject)
+        s["id"] = FIFF.FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE
+        srcs.append(s)
+
+    src = SourceSpaces(srcs, dict(working_dir=os.getcwd(), command_line="None"))
+    if add_dist:
+        add_source_space_distances(src, dist_limit=np.inf)
+    return src
+
+
 def _make_voxel_ras_trans(move, ras, voxel_size):
     """Make a transformation from MRI_VOXEL to MRI surface RAS (i.e. MRI)."""
     assert voxel_size.ndim == 1
@@ -2900,6 +3095,8 @@ def _get_hemi(s):
         return "lh", 0, s["id"]
     elif s["id"] == FIFF.FIFFV_MNE_SURF_RIGHT_HEMI:
         return "rh", 1, s["id"]
+    elif s["id"] >= FIFF.FIFFV_MNE_SURF_SUBCORTICAL_ID_BASE:
+        return s.get("seg_name", "subcortical"), None, s["id"]
     else:
         raise ValueError(f"unknown surface ID {s['id']}")
 

mne/surface.py

@@ -2100,6 +2100,23 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None, use_flying_edge
     return out
 
 
+def _keep_largest_component(rr, tris):
+    """Keep only the largest connected component of a triangulated mesh."""
+    from scipy.sparse.csgraph import connected_components
+
+    if len(tris) == 0:
+        return rr, tris
+    n_comp, labels = connected_components(mesh_edges(tris), directed=False)
+    if n_comp == 1:
+        return rr, tris
+    largest = np.argmax(np.bincount(labels))
+    keep = labels == largest
+    new_index = np.full(len(rr), -1, int)
+    new_index[keep] = np.arange(keep.sum())
+    tris = new_index[tris[keep[tris].all(axis=1)]]
+    return rr[keep], tris
+
+
 @verbose
 def _vtk_smooth(pd, smooth, *, verbose=None):
     _validate_type(smooth, "numeric", smooth)

tools/vulture_allowlist.py

@@ -57,6 +57,9 @@
 
 deep
 
+# not yet referenced elsewhere in the package
+setup_subcortical_source_space
+
 # Backward compat or rarely used
 RawFIF
 estimate_head_mri_t