ENH: Improve make_head_surface options

azure-pipelines.yml

@@ -110,7 +110,7 @@ stages:
     - bash: |
         set -e
         python -m pip install --progress-bar off --upgrade pip setuptools wheel codecov
-        python -m pip install --progress-bar off mne-qt-browser[opengl] vtk scikit-learn pytest-error-for-skips python-picard "PySide6!=6.3.0" qtpy
+        python -m pip install --progress-bar off mne-qt-browser[opengl] pyvista scikit-learn pytest-error-for-skips python-picard "PySide6!=6.3.0" qtpy
         python -m pip uninstall -yq mne
         python -m pip install --progress-bar off --upgrade -e .[test]
       displayName: 'Install dependencies with pip'

doc/changes/latest.inc

@@ -24,12 +24,12 @@ Current (1.1.dev0)
 Enhancements
 ~~~~~~~~~~~~
 
-- Add ``mne-icalabel`` to :func:`mne.sys_info` (:gh:`10615` by `Adam Li`_)
-
 - Add support for ahdr files in :func:`mne.io.read_raw_brainvision` (:gh:`10515` by :newcontrib:`Alessandro Tonin`)
 
 - Add support for reading data from Gowerlabs devices to :func:`mne.io.read_raw_snirf` (:gh:`10555` by :newcontrib:`Samuel Powell` and `Robert Luke`_)
 
+- Add ``mne-icalabel`` to :func:`mne.sys_info` (:gh:`10615` by `Adam Li`_)
+
 - Add support for ``overview_mode`` in :meth:`raw.plot() <mne.io.Raw.plot>` and related functions/methods (:gh:`10501` by `Eric Larson`_)
 
 - Add :meth:`mne.io.Raw.crop_by_annotations` method to get chunks of Raw data based on :class:`mne.Annotations`. (:gh:`10460` by `Alex Gramfort`_)
@@ -52,6 +52,8 @@ Enhancements
 
 - Add keyboard shortcuts to toggle :meth:`mne.preprocessing.ICA.plot_properties` topomap channel types ('t') and power spectral density log-scale ('l') (:gh:`10557` by `Alex Rockhill`_)
 
+- Add ``--mri``, and ``--threshold`` options to :ref:`mne make_scalp_surfaces` to improve head surface mesh extraction (:gh:`10591` by `Eric Larson`_)
+
 - Add :func:`mne.preprocessing.compute_bridged_electrodes` to detect EEG electrodes with shared spatial sources due to a conductive medium connecting two or more electrodes, add :ref:`ex-eeg-bridging` for an example and :func:`mne.viz.plot_bridged_electrodes` to help visualize (:gh:`10571` by `Alex Rockhill`_)
 
 - Add ``'voronoi'`` as an option for the ``image_interp`` argument in :func:`mne.viz.plot_topomap` to plot a topomap without interpolation using a Voronoi parcelation (:gh:`10571` by `Alex Rockhill`_)

doc/conf.py

@@ -308,7 +308,7 @@ def __call__(self, gallery_conf, fname, when):
         except ImportError:
             BackgroundPlotter = None  # noqa
         try:
-            from vtk import vtkPolyData  # noqa
+            from vtkmodules.vtkCommonDataModel import vtkPolyData  # noqa
         except ImportError:
             vtkPolyData = None  # noqa
         try:

mne/_freesurfer.py

@@ -683,9 +683,12 @@ def _get_head_surface(surf, subject, subjects_dir, bem=None, verbose=None):
         try_fnames = [op.join(subject_dir, 'bem', f'{subject}-head-dense.fif'),
                       op.join(subject_dir, 'surf', 'lh.seghead')]
     else:
-        try_fnames = [op.join(subject_dir, 'bem', 'outer_skin.surf'),
-                      op.join(subject_dir, 'bem', 'flash', 'outer_skin.surf'),
-                      op.join(subject_dir, 'bem', f'{subject}-head.fif')]
+        try_fnames = [
+            op.join(subject_dir, 'bem', 'outer_skin.surf'),
+            op.join(subject_dir, 'bem', 'flash', 'outer_skin.surf'),
+            op.join(subject_dir, 'bem', f'{subject}-head-sparse.fif'),
+            op.join(subject_dir, 'bem', f'{subject}-head.fif'),
+        ]
     for fname in try_fnames:
         if op.exists(fname):
             logger.info(f'Using {op.basename(fname)} for head surface.')

mne/bem.py

@@ -33,7 +33,7 @@
 from .utils import (verbose, logger, run_subprocess, get_subjects_dir, warn,
                     _pl, _validate_type, _TempDir, _check_freesurfer_home,
                     _check_fname, has_nibabel, _check_option, path_like,
-                    _on_missing, _import_h5io_funcs)
+                    _on_missing, _import_h5io_funcs, _ensure_int)
 
 
 # ############################################################################
@@ -261,10 +261,13 @@ def _check_complete_surface(surf, copy=False, incomplete='raise', extra=''):
     surf = complete_surface_info(surf, copy=copy, verbose=False)
     fewer = np.where([len(t) < 3 for t in surf['neighbor_tri']])[0]
     if len(fewer) > 0:
+        fewer = list(fewer)
+        fewer = (fewer[:80] + ['...']) if len(fewer) > 80 else fewer
+        fewer = ', '.join(str(f) for f in fewer)
         msg = ('Surface {} has topological defects: {:.0f} / {:.0f} vertices '
                'have fewer than three neighboring triangles [{}]{}'
-               .format(_bem_surf_name[surf['id']], len(fewer), surf['ntri'],
-                       ', '.join(str(f) for f in fewer), extra))
+               .format(_bem_surf_name[surf['id']], len(fewer), len(surf['rr']),
+                       fewer, extra))
         _on_missing(on_missing=incomplete, msg=msg, name='on_defects')
     return surf
 
@@ -480,7 +483,7 @@ def _surfaces_to_bem(surfs, ids, sigmas, ico=None, rescale=True,
                          'number of elements (1 or 3)')
     for si, surf in enumerate(surfs):
         if isinstance(surf, str):
-            surfs[si] = read_surface(surf, return_dict=True)[-1]
+            surfs[si] = surf = read_surface(surf, return_dict=True)[-1]
     # Downsampling if the surface is isomorphic with a subdivided icosahedron
     if ico is not None:
         for si, surf in enumerate(surfs):
@@ -1571,7 +1574,7 @@ def write_bem_surfaces(fname, surfs, overwrite=False, verbose=None):
 
 @verbose
 def write_head_bem(fname, rr, tris, on_defects='raise', overwrite=False,
-                   verbose=None):
+                   *, verbose=None):
     """Write a head surface to a fiff file.
 
     Parameters
@@ -2065,9 +2068,16 @@ def _check_file(fname, overwrite):
         raise IOError(f'File {fname} exists, use --overwrite to overwrite it')
 
 
+_tri_levels = dict(
+    medium=30000,
+    sparse=2500,
+)
+
+
 @verbose
 def make_scalp_surfaces(subject, subjects_dir=None, force=True,
-                        overwrite=False, no_decimate=False, verbose=None):
+                        overwrite=False, no_decimate=False, *,
+                        threshold=20, mri='T1.mgz', verbose=None):
     """Create surfaces of the scalp and neck.
 
     The scalp surfaces are required for using the MNE coregistration GUI, and
@@ -2080,29 +2090,35 @@ def make_scalp_surfaces(subject, subjects_dir=None, force=True,
     %(subjects_dir)s
     force : bool
         Force creation of the surface even if it has some topological defects.
-        Defaults to ``True``.
+        Defaults to ``True``. See :ref:`tut-fix-meshes` for ideas on how to
+        fix problematic meshes.
     %(overwrite)s
     no_decimate : bool
         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.
+    threshold : int
+        The threshold to use with the MRI in the call to ``mkheadsurf``.
+        The default is 20.
+
+        .. versionadded:: 1.1
+    mri : str
+        The MRI to use. Should exist in ``$SUBJECTS_DIR/$SUBJECT/mri``.
+
+        .. versionadded:: 1.1
     %(verbose)s
     """
-    this_env = deepcopy(os.environ)
     subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
-    this_env['SUBJECTS_DIR'] = subjects_dir
-    this_env['SUBJECT'] = subject
-    this_env['subjdir'] = subjects_dir + '/' + subject
-    if 'FREESURFER_HOME' not in this_env:
-        raise RuntimeError('The FreeSurfer environment needs to be set up '
-                           'for this script')
     incomplete = 'warn' if force else 'raise'
     subj_path = op.join(subjects_dir, subject)
     if not op.exists(subj_path):
         raise RuntimeError('%s does not exist. Please check your subject '
                            'directory path.' % subj_path)
 
-    mri = 'T1.mgz' if op.exists(op.join(subj_path, 'mri', 'T1.mgz')) else 'T1'
+    # Backward compat for old FreeSurfer (?)
+    _validate_type(mri, str, 'mri')
+    if mri == 'T1.mgz':
+        mri = mri if op.exists(op.join(subj_path, 'mri', mri)) else 'T1'
 
     logger.info('1. Creating a dense scalp tessellation with mkheadsurf...')
 
@@ -2116,9 +2132,21 @@ def check_seghead(surf_path=op.join(subj_path, 'surf')):
         return surf
 
     my_seghead = check_seghead()
+    threshold = _ensure_int(threshold, 'threshold')
     if my_seghead is None:
-        run_subprocess(['mkheadsurf', '-subjid', subject, '-srcvol', mri],
-                       env=this_env)
+        this_env = deepcopy(os.environ)
+        this_env['SUBJECTS_DIR'] = subjects_dir
+        this_env['SUBJECT'] = subject
+        this_env['subjdir'] = subjects_dir + '/' + subject
+        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)
 
     surf = check_seghead()
     if surf is None:
@@ -2133,18 +2161,20 @@ def check_seghead(surf_path=op.join(subj_path, 'surf')):
     logger.info('2. Creating %s ...' % dense_fname)
     _check_file(dense_fname, overwrite)
     # Helpful message if we get a topology error
-    msg = '\n\nConsider using --force as an additional input parameter.'
+    msg = ('\n\nConsider using pymeshfix directly to fix the mesh, or --force '
+           'to ignore the problem.')
     surf = _surfaces_to_bem(
         [surf], [FIFF.FIFFV_BEM_SURF_ID_HEAD], [1],
         incomplete=incomplete, extra=msg)[0]
     write_bem_surfaces(dense_fname, surf, overwrite=overwrite)
-    levels = 'medium', 'sparse'
-    tris = [] if no_decimate else [30000, 2500]
     if os.getenv('_MNE_TESTING_SCALP', 'false') == 'true':
         tris = [len(surf['tris'])]  # don't actually decimate
-    for ii, (n_tri, level) in enumerate(zip(tris, levels), 3):
-        logger.info('%i. Creating %s tessellation...' % (ii, level))
-        logger.info('%i.1 Decimating the dense tessellation...' % ii)
+    for ii, (level, n_tri) in enumerate(_tri_levels.items(), 3):
+        if no_decimate:
+            break
+        logger.info(f'{ii}. Creating {level} tessellation...')
+        logger.info(f'{ii}.1 Decimating the dense tessellation '
+                    f'({len(surf["tris"])} -> {n_tri} triangles)...')
         points, tris = decimate_surface(points=surf['rr'],
                                         triangles=surf['tris'],
                                         n_triangles=n_tri)
@@ -2156,3 +2186,4 @@ def check_seghead(surf_path=op.join(subj_path, 'surf')):
             [FIFF.FIFFV_BEM_SURF_ID_HEAD], [1], rescale=False,
             incomplete=incomplete, extra=msg)
         write_bem_surfaces(dec_fname, dec_surf, overwrite=overwrite)
+    logger.info('[done]')

mne/commands/mne_make_scalp_surfaces.py

@@ -34,9 +34,13 @@ def run():
                       help='Overwrite previously computed surface')
     parser.add_option('-s', '--subject', dest='subject',
                       help='The name of the subject', type='str')
+    parser.add_option('-m', '--mri', dest='mri', type='str', default='T1.mgz',
+                      help='The MRI file to process using mkheadsurf.')
     parser.add_option('-f', '--force', dest='force', action='store_true',
                       help='Force creation of the surface even if it has '
                            'some topological defects.')
+    parser.add_option('-t', '--threshold', dest='threshold', type='int',
+                      default=20, help='Threshold value to use with the MRI.')
     parser.add_option("-d", "--subjects-dir", dest="subjects_dir",
                       help="Subjects directory", default=subjects_dir)
     parser.add_option("-n", "--no-decimate", dest="no_decimate",
@@ -56,6 +60,8 @@ def run():
         force=options.force,
         overwrite=options.overwrite,
         no_decimate=options.no_decimate,
+        threshold=options.threshold,
+        mri=options.mri,
         verbose=options.verbose)
 
 

mne/commands/tests/test_commands.py

@@ -8,6 +8,7 @@
 import pytest
 from numpy.testing import assert_equal, assert_allclose
 
+import mne
 from mne import (concatenate_raws, read_bem_surfaces, read_surface,
                  read_source_spaces, read_bem_solution)
 from mne.bem import ConductorModel
@@ -22,7 +23,7 @@
                           mne_prepare_bem_model, mne_sys_info)
 from mne.datasets import testing
 from mne.io import read_raw_fif, read_info
-from mne.utils import (requires_mne, requires_vtk, requires_freesurfer,
+from mne.utils import (requires_mne, requires_freesurfer,
                        requires_nibabel, ArgvSetter,
                        _stamp_to_dt, _record_warnings)
 
@@ -133,11 +134,11 @@ def test_kit2fiff():
 
 @pytest.mark.slowtest
 @pytest.mark.ultraslowtest
-@requires_vtk
 @testing.requires_testing_data
 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
@@ -148,16 +149,18 @@ def test_make_scalp_surfaces(tmp_path, monkeypatch):
     os.mkdir(surf_path_new)
     subj_dir = op.join(tempdir, 'sample', 'bem')
     os.mkdir(subj_dir)
-    shutil.copy(op.join(surf_path, 'lh.seghead'), surf_path_new)
 
     cmd = ('-s', 'sample', '--subjects-dir', tempdir)
-    monkeypatch.setenv('_MNE_TESTING_SCALP', 'true')
+    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')
     with ArgvSetter(cmd, disable_stdout=False, disable_stderr=False):
         monkeypatch.delenv('FREESURFER_HOME', None)
         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(dense_fname)

mne/coreg.py

@@ -53,10 +53,12 @@
 surf_dirname = os.path.join(subject_dirname, 'surf')
 bem_fname = os.path.join(bem_dirname, "{subject}-{name}.fif")
 head_bem_fname = pformat(bem_fname, name='head')
+head_sparse_fname = pformat(bem_fname, name='head-sparse')
 fid_fname = pformat(bem_fname, name='fiducials')
 fid_fname_general = os.path.join(bem_dirname, "{head}-fiducials.fif")
 src_fname = os.path.join(bem_dirname, '{subject}-{spacing}-src.fif')
 _head_fnames = (os.path.join(bem_dirname, 'outer_skin.surf'),
+                head_sparse_fname,
                 head_bem_fname)
 _high_res_head_fnames = (os.path.join(bem_dirname, '{subject}-head-dense.fif'),
                          os.path.join(surf_dirname, 'lh.seghead'),