Fixes make_scalp_surfaces (#13024)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Marijn van Vliet <w.m.vanvliet@gmail.com>
Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

Author: 4 people

doc/changes/dev/13024.bugfix.rst

@@ -0,0 +1 @@
+Fix bug where outdated ``seghead.mgz`` files were reused in :func:`mne.bem.make_scalp_surfaces`. Add a new parameter ``reuse_seghead`` to control whether to reuse existing ``seghead.mgz`` files, by `Victor Ferat`_.

mne/bem.py

@@ -2349,6 +2349,7 @@ def make_scalp_surfaces(
     force=True,
     overwrite=False,
     no_decimate=False,
+    reuse_seghead=False,
     *,
     threshold=20,
     mri="T1.mgz",
@@ -2373,6 +2374,12 @@ def make_scalp_surfaces(
         Disable the "medium" and "sparse" decimations. In this case, only
         a "dense" surface will be generated. Defaults to ``False``, i.e.,
         create surfaces for all three types of decimations.
+    reuse_seghead : bool
+        Whether to reuse existing head segmentation files. If ``True``,
+        the existing files will be used if they exist. If ``False``
+        (default), the head segmentation will be recomputed.
+
+        .. versionadded:: 1.12
     threshold : int
         The threshold to use with the MRI in the call to ``mkheadsurf``.
         The default is ``20``.
@@ -2397,30 +2404,58 @@ def make_scalp_surfaces(
     if mri == "T1.mgz":
         mri = mri if (subj_path / "mri" / mri).exists() else "T1"
 
-    logger.info("1. Creating a dense scalp tessellation with mkheadsurf...")
+    threshold = _ensure_int(threshold, "threshold")
 
-    def check_seghead(surf_path=subj_path / "surf"):
-        surf = None
-        for k in ["lh.seghead", "lh.smseghead"]:
-            this_surf = surf_path / k
-            if this_surf.exists():
-                surf = this_surf
-                break
-        return surf
+    # Check for existing files
+    seghead_mgz_path = subj_path / "mri" / "seghead.mgz"
+    seghead_surf_path = subj_path / "surf" / "lh.seghead"
+    smseghead_surf_path = subj_path / "surf" / "lh.smseghead"
 
-    my_seghead = check_seghead()
-    threshold = _ensure_int(threshold, "threshold")
-    if my_seghead is None:
+    bem_dir = subjects_dir / subject / "bem"
+    fname_template = bem_dir / (f"{subject}-head-{{}}.fif")
+    dense_fname = str(fname_template).format("dense")
+    _check_file(dense_fname, overwrite)
+
+    for level in _tri_levels:
+        dec_fname = str(fname_template).format(level)
+        if overwrite:
+            if os.path.exists(dec_fname):
+                logger.info(f"Removing previously existing {dec_fname}.")
+                os.remove(dec_fname)
+        else:
+            if no_decimate:
+                if os.path.exists(dec_fname):
+                    raise OSError(
+                        f"Trying to generate new scalp surfaces"
+                        f"but {dec_fname} already exists."
+                        f"To avoid mixing different scalp surface solutions, "
+                        f"delete this file or use overwrite to automatically delete it."
+                    )
+            else:
+                _check_file(dec_fname, overwrite)
+
+    _check_freesurfer_home()
+
+    if reuse_seghead:
+        if seghead_surf_path.exists():
+            surf = seghead_surf_path
+        elif smseghead_surf_path.exists():
+            surf = smseghead_surf_path
+        else:
+            raise ValueError(
+                "No existing scalp surface found. Please check your subject's surf "
+                "folder or set reuse_seghead to False to recompute the surfaces."
+            )
+        logger.info(f"1. Using existing scalp tessellation {surf} ...")
+    else:
+        _check_file(seghead_mgz_path, overwrite)
+        _check_file(seghead_surf_path, overwrite)
+        _check_file(smseghead_surf_path, overwrite)
+        logger.info("1. Creating a dense scalp tessellation with mkheadsurf...")
         this_env = deepcopy(os.environ)
         this_env["SUBJECTS_DIR"] = str(subjects_dir)
         this_env["SUBJECT"] = subject
         this_env["subjdir"] = str(subj_path)
-        if "FREESURFER_HOME" not in this_env:
-            raise RuntimeError(
-                "The FreeSurfer environment needs to be set up to use "
-                "make_scalp_surfaces to create the outer skin surface "
-                "lh.seghead"
-            )
         run_subprocess(
             [
                 "mkheadsurf",
@@ -2435,18 +2470,16 @@ def check_seghead(surf_path=subj_path / "surf"):
             ],
             env=this_env,
         )
+        if os.path.exists(seghead_surf_path):
+            surf = seghead_surf_path
+        elif os.path.exists(smseghead_surf_path):
+            surf = smseghead_surf_path
+        else:
+            raise ValueError("mkheadsurf did not produce the standard output file.")
 
-    surf = check_seghead()
-    if surf is None:
-        raise RuntimeError("mkheadsurf did not produce the standard output file.")
-
-    bem_dir = subjects_dir / subject / "bem"
-    if not bem_dir.is_dir():
-        os.mkdir(bem_dir)
-    fname_template = bem_dir / (f"{subject}-head-{{}}.fif")
-    dense_fname = str(fname_template).format("dense")
     logger.info(f"2. Creating {dense_fname} ...")
-    _check_file(dense_fname, overwrite)
+    bem_dir.mkdir(exist_ok=True)
+
     # Helpful message if we get a topology error
     msg = (
         "\n\nConsider using pymeshfix directly to fix the mesh, or --force "
@@ -2471,7 +2504,6 @@ def check_seghead(surf_path=subj_path / "surf"):
         )
         dec_fname = str(fname_template).format(level)
         logger.info(f"{ii}.2 Creating {dec_fname}")
-        _check_file(dec_fname, overwrite)
         dec_surf = _surfaces_to_bem(
             [dict(rr=points, tris=tris)],
             [FIFF.FIFFV_BEM_SURF_ID_HEAD],

mne/commands/mne_make_scalp_surfaces.py

@@ -73,6 +73,13 @@ def run():
         help="Disable medium and sparse decimations (dense only)",
         action="store_true",
     )
+    parser.add_option(
+        "-r",
+        "--reuse-seghead",
+        dest="reuse_seghead",
+        action="store_true",
+        help="Whether to reuse existing head segmentation files.",
+    )
     _add_verbose_flag(parser)
     options, args = parser.parse_args()
 
@@ -89,6 +96,7 @@ def run():
         no_decimate=options.no_decimate,
         threshold=options.threshold,
         mri=options.mri,
+        reuse_seghead=options.reuse_seghead,
         verbose=options.verbose,
     )
 

mne/commands/tests/test_commands.py

@@ -170,51 +170,129 @@ def test_kit2fiff():
     check_usage(mne_kit2fiff, force_help=True)
 
 
-@pytest.mark.slowtest
 @pytest.mark.ultraslowtest
 @testing.requires_testing_data
+@requires_freesurfer("mkheadsurf")
 def test_make_scalp_surfaces(tmp_path, monkeypatch):
     """Test mne make_scalp_surfaces."""
     pytest.importorskip("nibabel")
     pytest.importorskip("pyvista")
     check_usage(mne_make_scalp_surfaces)
     has = "SUBJECTS_DIR" in os.environ
-    # Copy necessary files to avoid FreeSurfer call
-    tempdir = str(tmp_path)
-    surf_path = op.join(subjects_dir, "sample", "surf")
-    surf_path_new = op.join(tempdir, "sample", "surf")
-    os.mkdir(op.join(tempdir, "sample"))
-    os.mkdir(surf_path_new)
-    subj_dir = op.join(tempdir, "sample", "bem")
-    os.mkdir(subj_dir)
+    freesurfer_home = os.environ.get("FREESURFER_HOME")
 
-    cmd = ("-s", "sample", "--subjects-dir", tempdir)
     monkeypatch.setattr(
         mne.bem,
         "decimate_surface",
         lambda points, triangles, n_triangles: (points, triangles),
     )
-    dense_fname = op.join(subj_dir, "sample-head-dense.fif")
-    medium_fname = op.join(subj_dir, "sample-head-medium.fif")
+
+    tempdir = str(tmp_path)
+    t1_path = op.join(subjects_dir, "sample", "mri", "T1.mgz")
+    t1_path_new = op.join(tempdir, "sample", "mri", "T1.mgz")
+
+    headseg_path = op.join(subjects_dir, "sample", "mri", "seghead.mgz")
+    headseg_path_new = op.join(tempdir, "sample", "mri", "seghead.mgz")
+
+    surf_path = op.join(subjects_dir, "sample", "surf", "lh.seghead")
+    surf_path_new = op.join(tempdir, "sample", "surf", "lh.seghead")
+
+    dense_fname = op.join(tempdir, "sample", "bem", "sample-head-dense.fif")
+    medium_fname = op.join(tempdir, "sample", "bem", "sample-head-medium.fif")
+    sparse_fname = op.join(tempdir, "sample", "bem", "sample-head-sparse.fif")
+
+    os.makedirs(op.join(tempdir, "sample", "mri"), exist_ok=True)
+    os.makedirs(op.join(tempdir, "sample", "surf"), exist_ok=True)
+
+    shutil.copy(t1_path, t1_path_new)
+    shutil.copy(headseg_path, headseg_path_new)
+    shutil.copy(surf_path, surf_path_new)
+
+    cmd = (
+        "-s",
+        "sample",
+        "--subjects-dir",
+        tempdir,
+        "--no-decimate",
+        "--reuse-seghead",
+    )
     with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
         monkeypatch.delenv("FREESURFER_HOME", raising=False)
-        with pytest.raises(RuntimeError, match="The FreeSurfer environ"):
+        with pytest.raises(RuntimeError, match="The FREESURFER_HOME environment"):
             mne_make_scalp_surfaces.run()
-        shutil.copy(op.join(surf_path, "lh.seghead"), surf_path_new)
-        monkeypatch.setenv("FREESURFER_HOME", tempdir)
+
+        monkeypatch.setenv("FREESURFER_HOME", freesurfer_home)
         mne_make_scalp_surfaces.run()
+
+        assert op.isfile(headseg_path_new)
+        assert op.isfile(surf_path_new)
+        assert op.isfile(dense_fname)
+        assert not op.isfile(medium_fname)
+        assert not op.isfile(sparse_fname)
+
+        # actually check the outputs
+        head_py = read_bem_surfaces(dense_fname)
+        assert_equal(len(head_py), 1)
+        head_py = head_py[0]
+        head_c = read_bem_surfaces(
+            op.join(subjects_dir, "sample", "bem", "sample-head-dense.fif")
+        )[0]
+        assert_allclose(head_py["rr"], head_c["rr"])
+
+    cmd = ("-s", "sample", "--subjects-dir", tempdir)
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        with pytest.raises(OSError, match="use --overwrite to overwrite it"):
+            mne_make_scalp_surfaces.run()
+
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--overwrite")
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        mne_make_scalp_surfaces.run()
+        assert op.isfile(headseg_path_new)
+        assert op.isfile(surf_path_new)
         assert op.isfile(dense_fname)
         assert op.isfile(medium_fname)
-        with pytest.raises(OSError, match="overwrite"):
+        assert op.isfile(sparse_fname)
+
+    os.remove(headseg_path_new)
+    os.remove(surf_path_new)
+    os.remove(dense_fname)
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--no-decimate")
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        with pytest.raises(OSError, match="Trying to generate new scalp surfaces"):
+            mne_make_scalp_surfaces.run()
+
+    os.remove(medium_fname)
+    os.remove(sparse_fname)
+    cmd = (
+        "-s",
+        "sample",
+        "--subjects-dir",
+        tempdir,
+        "--no-decimate",
+        "--reuse-seghead",
+    )
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        with pytest.raises(ValueError, match="No existing scalp surface found"):
             mne_make_scalp_surfaces.run()
-    # actually check the outputs
-    head_py = read_bem_surfaces(dense_fname)
-    assert_equal(len(head_py), 1)
-    head_py = head_py[0]
-    head_c = read_bem_surfaces(
-        op.join(subjects_dir, "sample", "bem", "sample-head-dense.fif")
-    )[0]
-    assert_allclose(head_py["rr"], head_c["rr"])
+
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--no-decimate", "--overwrite")
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        mne_make_scalp_surfaces.run()
+        assert op.isfile(headseg_path_new)
+        assert op.isfile(surf_path_new)
+        assert op.isfile(dense_fname)
+        assert not op.isfile(medium_fname)
+        assert not op.isfile(sparse_fname)
+
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--no-decimate", "--overwrite")
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        mne_make_scalp_surfaces.run()
+        assert op.isfile(headseg_path_new)
+        assert op.isfile(surf_path_new)
+        assert op.isfile(dense_fname)
+        assert not op.isfile(medium_fname)
+        assert not op.isfile(sparse_fname)
+
     if not has:
         assert "SUBJECTS_DIR" not in os.environ
 

mne/tests/test_bem.py

@@ -43,7 +43,12 @@
 from mne.io import read_info
 from mne.surface import _get_ico_surface, read_surface
 from mne.transforms import translation
-from mne.utils import _record_warnings, catch_logging, check_version
+from mne.utils import (
+    _record_warnings,
+    catch_logging,
+    check_version,
+    requires_freesurfer,
+)
 
 fname_raw = Path(__file__).parents[1] / "io" / "tests" / "data" / "test_raw.fif"
 subjects_dir = testing.data_path(download=False) / "subjects"
@@ -518,11 +523,27 @@ def test_io_head_bem(tmp_path):
     assert np.allclose(head["tris"], head_defect["tris"])
 
 
-@pytest.mark.slowtest  # ~4 s locally
+@pytest.mark.slowtest
+@requires_freesurfer("mkheadsurf")
+@testing.requires_testing_data
 def test_make_scalp_surfaces_topology(tmp_path, monkeypatch):
-    """Test topology checks for make_scalp_surfaces."""
+    """Test make_scalp_surfaces and topology checks."""
     pytest.importorskip("pyvista")
     pytest.importorskip("nibabel")
+
+    # tests on 'sample'
+    subject = "sample"
+    subjects_dir = testing.data_path(download=False)
+    with pytest.raises(OSError, match="use --overwrite to overwrite it"):
+        make_scalp_surfaces(
+            subject, subjects_dir, force=False, verbose=True, overwrite=False
+        )
+
+    make_scalp_surfaces(
+        subject, subjects_dir, force=False, verbose=True, overwrite=True
+    )
+
+    # tests on custom surface
     subjects_dir = tmp_path
     subject = "test"
     surf_dir = subjects_dir / subject / "surf"
@@ -550,6 +571,7 @@ def _decimate_surface(points, triangles, n_triangles):
         make_scalp_surfaces(
             subject, subjects_dir, force=False, verbose=True, overwrite=True
         )
+
     bem_dir = subjects_dir / subject / "bem"
     sparse_path = bem_dir / f"{subject}-head-sparse.fif"
     assert not sparse_path.is_file()