Fixes make_scalp_surfaces

mne/bem.py

@@ -2393,56 +2393,75 @@ 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
+    _check_file(subj_path / "mri" / "seghead.mgz", overwrite)
 
-    my_seghead = check_seghead()
-    threshold = _ensure_int(threshold, "threshold")
-    if my_seghead is None:
-        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",
-                "-subjid",
-                subject,
-                "-srcvol",
-                mri,
-                "-thresh1",
-                str(threshold),
-                "-thresh2",
-                str(threshold),
-            ],
-            env=this_env,
-        )
+    seghead_path = subj_path / "surf" / "lh.seghead"
+    _check_file(seghead_path, overwrite)
 
-    surf = check_seghead()
-    if surf is None:
-        raise RuntimeError("mkheadsurf did not produce the standard output file.")
+    smseghead_path = subj_path / "surf" / "lh.smseghead"
+    _check_file(smseghead_path, overwrite)
 
     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)
+
+    for level in _tri_levels.items():
+        dec_fname = str(fname_template).format(level)
+        if overwrite:
+            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)
+
+    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",
+            "-subjid",
+            subject,
+            "-srcvol",
+            mri,
+            "-thresh1",
+            str(threshold),
+            "-thresh2",
+            str(threshold),
+        ],
+        env=this_env,
+    )
+    if os.path.exists(seghead_path):
+        surf = seghead_path
+    elif os.path.exists(smseghead_path):
+        surf = smseghead_path
+    else:
+        raise ValueError("mkheadsurf did not produce the standard output file.")
+
+    logger.info(f"2. Creating {dense_fname} ...")
+    if not bem_dir.is_dir():
+        os.mkdir(bem_dir)
+
     # Helpful message if we get a topology error
     msg = (
         "\n\nConsider using pymeshfix directly to fix the mesh, or --force "
@@ -2452,6 +2471,7 @@ def check_seghead(surf_path=subj_path / "surf"):
         [surf], [FIFF.FIFFV_BEM_SURF_ID_HEAD], [1], incomplete=incomplete, extra=msg
     )[0]
     write_bem_surfaces(dense_fname, surf, overwrite=overwrite)
+
     if os.getenv("_MNE_TESTING_SCALP", "false") == "true":
         tris = [len(surf["tris"])]  # don't actually decimate
     for ii, (level, n_tri) in enumerate(_tri_levels.items(), 3):
@@ -2467,7 +2487,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/tests/test_commands.py

@@ -12,7 +12,6 @@
 import pytest
 from numpy.testing import assert_allclose, assert_equal
 
-import mne
 from mne import (
     concatenate_raws,
     read_bem_solution,
@@ -167,53 +166,74 @@ 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)
+    t1_path = op.join(subjects_dir, "sample", "mri", "T1.mgz")
+    t1_path_new = op.join(tempdir, "sample", "mri", "T1.mgz")
 
-    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")
+    headseg_path = op.join(tempdir, "sample", "mri", "seghead.mgz")
+    surf_path = 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)
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--no-decimate")
     with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
         monkeypatch.delenv("FREESURFER_HOME", raising=False)
         with pytest.raises(RuntimeError, match="The FreeSurfer environ"):
             mne_make_scalp_surfaces.run()
-        shutil.copy(op.join(surf_path, "lh.seghead"), surf_path_new)
+
         monkeypatch.setenv("FREESURFER_HOME", tempdir)
         mne_make_scalp_surfaces.run()
+
+        assert op.isfile(headseg_path)
+        assert op.isfile(surf_path)
+        assert op.isfile(dense_fname)
+        assert not op.isfile(medium_fname)
+        assert not op.isfile(sparse_fname)
+
+    cmd = ("-s", "sample", "--subjects-dir", tempdir)
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        with pytest.raises(RuntimeError, 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)
+        assert op.isfile(surf_path)
         assert op.isfile(dense_fname)
         assert op.isfile(medium_fname)
-        with pytest.raises(OSError, match="overwrite"):
+        assert 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)
+        assert op.isfile(surf_path)
+        assert op.isfile(dense_fname)
+        assert not op.isfile(medium_fname)
+        assert not op.isfile(sparse_fname)
+
+    os.remove(headseg_path)
+    os.remove(surf_path)
+    os.remove(dense_fname)
+    cmd = ("-s", "sample", "--subjects-dir", tempdir, "--no-decimate")
+    with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
+        with pytest.raises(RuntimeError, match="Trying to generate new scalp surfaces"):
             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"])
-    if not has:
-        assert "SUBJECTS_DIR" not in os.environ
 
 
 @pytest.mark.slowtest

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"
@@ -519,11 +524,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"
@@ -551,6 +572,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()