diff --git a/doc/api/forward.rst b/doc/api/forward.rst index 5abcd5178fc..e6cb9c681d9 100644 --- a/doc/api/forward.rst +++ b/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 diff --git a/mne/__init__.pyi b/mne/__init__.pyi index 9c5573e0686..eca3d5f3e9c 100644 --- a/mne/__init__.pyi +++ b/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, ) diff --git a/mne/_fiff/constants.py b/mne/_fiff/constants.py index aced3454d57..5ebf2f0b3e0 100644 --- a/mne/_fiff/constants.py +++ b/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)) # diff --git a/mne/_fiff/tests/test_constants.py b/mne/_fiff/tests/test_constants.py index 102510ab33c..eea628ed38c 100644 --- a/mne/_fiff/tests/test_constants.py +++ b/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", diff --git a/mne/source_estimate.py b/mne/source_estimate.py index edb1ae45764..ccabccd14e3 100644 --- a/mne/source_estimate.py +++ b/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 diff --git a/mne/source_space/__init__.pyi b/mne/source_space/__init__.pyi index aeb7657bd33..9031e2bc2cc 100644 --- a/mne/source_space/__init__.pyi +++ b/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, ) diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 3bf6ec92bbf..6f6cf894369 100644 --- a/mne/source_space/_source_space.py +++ b/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']}") diff --git a/mne/surface.py b/mne/surface.py index c3075bf7d9d..b163fe8ca50 100644 --- a/mne/surface.py +++ b/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) diff --git a/tools/vulture_allowlist.py b/tools/vulture_allowlist.py index 1367a2cac25..2f890b2cdfe 100644 --- a/tools/vulture_allowlist.py +++ b/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