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/doc/changes/dev/14130.newfeature.rst b/doc/changes/dev/14130.newfeature.rst new file mode 100644 index 00000000000..01b5859d65a --- /dev/null +++ b/doc/changes/dev/14130.newfeature.rst @@ -0,0 +1 @@ +Add :func:`mne.setup_subcortical_source_space` to build a :class:`~mne.SourceSpaces` from a triangulated surface mesh of a subcortical or cerebellar structure, either from a ``label`` or externally-produced mesh ``surface``, by `Payam Sadeghi-Shabestari`_. diff --git a/mne/__init__.pyi b/mne/__init__.pyi index 1b5821536a5..0d09fb17098 100644 --- a/mne/__init__.pyi +++ b/mne/__init__.pyi @@ -165,6 +165,7 @@ __all__ = [ "set_log_level", "set_memmap_min_size", "setup_source_space", + "setup_subcortical_source_space", "setup_volume_source_space", "simulation", "source_space", @@ -410,6 +411,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..89150756542 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_OFFSET = 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..fbf36bd8ac6 100644 --- a/mne/_fiff/tests/test_constants.py +++ b/mne/_fiff/tests/test_constants.py @@ -28,7 +28,7 @@ # https://github.com/mne-tools/fiff-constants/commits/master REPO = "mne-tools" -COMMIT = "9ccb09d69daa8332f2e7252638ba397b60ba2502" +COMMIT = "c434349e2961df29d0938e5ec8b522d0dca4efa0" # These are oddities that we won't address: iod_dups = (355, 359) # these are in both MEGIN and MNE files @@ -57,7 +57,7 @@ "viewkeys", "viewvalues", # Py2 ) -_tag_ignore_names = () # for fiff-constants pending updates +_tag_ignore_names = () _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 58af13de0c3..6ba1e9ae832 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -408,7 +408,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 @@ -436,7 +436,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 4d091988b41..c5c04b94c82 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -32,6 +32,7 @@ ) from .._freesurfer import ( _check_mri, + _get_aseg, _get_atlas_values, _get_mri_info_data, get_volume_labels_from_aseg, @@ -47,6 +48,8 @@ _create_surf_spacing, _get_ico_surface, _get_surf_neighbors, + _keep_largest_component, + _marching_cubes, _normalize_vectors, _tessellate_sphere_surf, _triangle_neighbors, @@ -295,6 +298,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)) @@ -303,19 +307,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_OFFSET + 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 @@ -1065,6 +1090,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 @@ -1461,7 +1487,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"]) @@ -2011,6 +2037,182 @@ 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. + + .. warning:: + This is **experimental** functionality. :class:`~mne.SourceSpaces` + created by this function are not (yet) compatible with morphing + (:class:`~mne.SourceMorph`), :func:`mne.extract_label_time_course` + does not yet know how to select vertices within a subcortical-surface + label, and plotting support is limited (for example, the + :meth:`~mne.MixedSourceEstimate.plot` method does not yet support + these source spaces). + + 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. See the warning above + for known gaps, to be addressed in follow-up work. + + .. 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_OFFSET + 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_OFFSET + 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 @@ -2948,6 +3150,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_OFFSET: + return s.get("seg_name", "subcortical"), None, s["id"] else: raise ValueError(f"unknown surface ID {s['id']}") diff --git a/mne/source_space/tests/test_source_space.py b/mne/source_space/tests/test_source_space.py index 7dbde2a7034..3fb7e06d3ab 100644 --- a/mne/source_space/tests/test_source_space.py +++ b/mne/source_space/tests/test_source_space.py @@ -26,8 +26,10 @@ read_bem_surfaces, read_freesurfer_lut, read_source_spaces, + read_surface, read_trans, setup_source_space, + setup_subcortical_source_space, setup_volume_source_space, write_source_spaces, ) @@ -674,6 +676,71 @@ def test_source_space_from_label(tmp_path, pass_ids): _compare_source_spaces(src, src_from_file, mode="approx") +@testing.requires_testing_data +def test_setup_subcortical_source_space(tmp_path): + """Test setting up a subcortical/cerebellar surface source space.""" + pytest.importorskip("nibabel") + pytest.importorskip("pyvista") + atlas_ids, _ = read_freesurfer_lut() + fname_surf = subjects_dir / "sample" / "surf" / "lh.white" + + # exactly one of label/surface must be given + with pytest.raises(ValueError, match="Exactly one of"): + setup_subcortical_source_space("sample", subjects_dir=subjects_dir) + with pytest.raises(ValueError, match="Exactly one of"): + setup_subcortical_source_space( + "sample", + label="Left-Amygdala", + surface=fname_surf, + subjects_dir=subjects_dir, + ) + + # label input: tessellate regions from the anatomical segmentation + labels = ["Left-Amygdala", "Left-Hippocampus"] + src_label = setup_subcortical_source_space( + "sample", label=labels, aseg="aseg", subjects_dir=subjects_dir + ) + assert src_label.kind == "subcortical_surf" + assert len(src_label) == len(labels) + for s, name in zip(src_label, labels): + assert s["type"] == "surf" + assert s["seg_name"] == name + assert s["id"] == FIFF.FIFFV_MNE_SURF_SUBCORTICAL_OFFSET + atlas_ids[name] + assert s["rr"].shape[1] == 3 and s["ntri"] > 0 + + # mesh input: an externally-supplied surface (path or rr/tris dict) + rr, tris = read_surface(fname_surf)[:2] + src_surface = setup_subcortical_source_space( + "sample", + surface=fname_surf, + subjects_dir=subjects_dir, + keep_largest_component=False, + ) + assert len(src_surface) == 1 + assert src_surface[0]["id"] == FIFF.FIFFV_MNE_SURF_SUBCORTICAL_OFFSET + assert src_surface[0]["rr"].shape[0] == rr.shape[0] + assert src_surface[0]["ntri"] == len(tris) + + src_surface_dict = setup_subcortical_source_space( + "sample", + surface=dict(rr=rr, tris=tris), + subjects_dir=subjects_dir, + keep_largest_component=False, + ) + assert_allclose(src_surface_dict[0]["rr"], src_surface[0]["rr"]) + + # I/O roundtrip + fname_temp = tmp_path / "subcortical-src.fif" + write_source_spaces(fname_temp, src_label) + src_read = read_source_spaces(fname_temp) + assert len(src_read) == len(src_label) + for orig, read in zip(src_label, src_read): + assert orig["seg_name"] == read["seg_name"] + assert orig["id"] == read["id"] + assert_allclose(orig["rr"], read["rr"], atol=1e-6) + assert_array_equal(orig["tris"], read["tris"]) + + @pytest.mark.slowtest @testing.requires_testing_data def test_source_space_exclusive_complete(src_volume_labels): diff --git a/mne/surface.py b/mne/surface.py index 52b6cb754c9..e284632b450 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1929,6 +1929,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 6625c2494f5..d0ee27b8177 100644 --- a/tools/vulture_allowlist.py +++ b/tools/vulture_allowlist.py @@ -62,6 +62,9 @@ deep +# not yet referenced elsewhere in the package +setup_subcortical_source_space + # Module-level __getattr__ (PEP 562), used by mne/surface.py and # mne/transforms.py to re-export their numba helpers lazily __getattr__