diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 915b05e31af..d7cf4232db4 100644 --- a/azure-pipelines.yml +++ b/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' diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 0f37bc13693..472fc247f3d 100644 --- a/doc/changes/latest.inc +++ b/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() ` 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`_) diff --git a/doc/conf.py b/doc/conf.py index 918836f4d3b..7f4bcc67cd5 100644 --- a/doc/conf.py +++ b/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: diff --git a/mne/_freesurfer.py b/mne/_freesurfer.py index 7b5e33f808d..8e643e96db5 100644 --- a/mne/_freesurfer.py +++ b/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.') diff --git a/mne/bem.py b/mne/bem.py index e419302fd3c..e53447933b0 100644 --- a/mne/bem.py +++ b/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]') diff --git a/mne/commands/mne_make_scalp_surfaces.py b/mne/commands/mne_make_scalp_surfaces.py index 1e0e4cbcd8a..9da7941384c 100644 --- a/mne/commands/mne_make_scalp_surfaces.py +++ b/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) diff --git a/mne/commands/tests/test_commands.py b/mne/commands/tests/test_commands.py index ed8de974470..2eb855d71c9 100644 --- a/mne/commands/tests/test_commands.py +++ b/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) diff --git a/mne/coreg.py b/mne/coreg.py index 1fcd814e6da..897fd5e9d66 100644 --- a/mne/coreg.py +++ b/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'), diff --git a/mne/surface.py b/mne/surface.py index 5bfd68fde89..0c395a6b6b0 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -572,17 +572,27 @@ def _points_outside_surface(rr, surf, n_jobs=None, verbose=None): return np.abs(np.sum(tot_angles, axis=0) / (2 * np.pi) - 1.0) > 1e-5 -def _surface_to_polydata(rr, tris=None): +def _surface_to_polydata(surf): import pyvista as pv - vertices = np.array(rr) - if tris is None: + vertices = np.array(surf['rr']) + if 'tris' not in surf: return pv.PolyData(vertices) else: - triangles = np.array(tris) + triangles = np.array(surf['tris']) triangles = np.c_[np.full(len(triangles), 3), triangles] return pv.PolyData(vertices, triangles) +def _polydata_to_surface(pd, normals=True): + from pyvista import PolyData + if not isinstance(pd, PolyData): + pd = PolyData(pd) + out = dict(rr=pd.points, tris=pd.faces.reshape(-1, 4)[:, 1:]) + if normals: + out['nn'] = pd.point_normals + return out + + class _CheckInside(object): """Efficiently check if points are inside a surface.""" @@ -617,13 +627,9 @@ def _init_old(self): def _init_pyvista(self): if not isinstance(self.surf, dict): self.pdata = self.surf - self.surf = dict( - rr=self.pdata.points, - tris=self.pdata.faces.reshape(-1, 4)[:, 1:], - nn=self.pdata.point_normals) + self.surf = _polydata_to_surface(self.pdata) else: - self.pdata = _surface_to_polydata( - self.surf['rr'], self.surf['tris']).clean() + self.pdata = _surface_to_polydata(self.surf).clean() @verbose def __call__(self, rr, n_jobs=None, verbose=None): @@ -643,7 +649,7 @@ def __call__(self, rr, n_jobs=None, verbose=None): return inside def _call_pyvista(self, rr): - pdata = _surface_to_polydata(rr) + pdata = _surface_to_polydata(dict(rr=rr)) out = pdata.select_enclosed_points(self.pdata, check_surface=False) return out['SelectedPoints'].astype(bool) @@ -818,11 +824,14 @@ def read_surface(fname, read_metadata=False, return_dict=False, ret += (dict(),) if return_dict: - ret += (dict(rr=ret[0], tris=ret[1], ntri=len(ret[1]), use_tris=ret[1], - np=len(ret[0])),) + ret += (_rr_tris_dict(ret[0], ret[1]),) return ret +def _rr_tris_dict(rr, tris): + return dict(rr=rr, tris=tris, ntri=len(tris), use_tris=tris, np=len(rr)) + + def _read_mri_surface(fname): surf = read_surface(fname, return_dict=True)[2] surf['rr'] /= 1000. @@ -1228,11 +1237,11 @@ def write_surface(fname, coords, faces, create_stamp='', volume_info=None, def _decimate_surface_vtk(points, triangles, n_triangles): """Aux function.""" try: - from vtk.util.numpy_support import \ + from vtkmodules.util.numpy_support import \ numpy_to_vtk, numpy_to_vtkIdTypeArray - from vtk.numpy_interface.dataset_adapter import WrapDataObject - from vtk import \ - vtkPolyData, vtkQuadricDecimation, vtkPoints, vtkCellArray + from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkCellArray + from vtkmodules.vtkCommonCore import vtkPoints + from vtkmodules.vtkFiltersCore import vtkQuadricDecimation except ImportError: raise ValueError('This function requires the VTK package to be ' 'installed') @@ -1261,10 +1270,9 @@ def _decimate_surface_vtk(points, triangles, n_triangles): reduction = 1 - (float(n_triangles) / len(triangles)) decimate.SetTargetReduction(reduction) decimate.Update() - out = WrapDataObject(decimate.GetOutput()) - rrs = out.Points - tris = out.Polygons.reshape(-1, 4)[:, 1:] - return rrs, tris + + out = _polydata_to_surface(decimate.GetOutput(), normals=False) + return out['rr'], out['tris'] def _decimate_surface_sphere(rr, tris, n_triangles): @@ -1317,7 +1325,7 @@ def _decimate_surface_sphere(rr, tris, n_triangles): @verbose -def decimate_surface(points, triangles, n_triangles, method='quadric', +def decimate_surface(points, triangles, n_triangles, method='quadric', *, verbose=None): """Decimate surface data. @@ -1754,16 +1762,13 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None): # Also vtkDiscreteFlyingEdges3D should be faster. # If we ever want not-discrete (continuous/float) marching cubes, # we should probably use vtkFlyingEdges3D rather than vtkMarchingCubes. - from vtk import (vtkImageData, vtkThreshold, - vtkWindowedSincPolyDataFilter, vtkDiscreteFlyingEdges3D, - vtkGeometryFilter, vtkDataSetAttributes, VTK_DOUBLE) - from vtk.util import numpy_support + from vtkmodules.vtkCommonDataModel import \ + vtkImageData, vtkDataSetAttributes + from vtkmodules.vtkFiltersCore import vtkThreshold + from vtkmodules.vtkFiltersGeneral import vtkDiscreteFlyingEdges3D + from vtkmodules.vtkFiltersGeometry import vtkGeometryFilter + from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk from scipy.ndimage import binary_dilation - _validate_type(smooth, 'numeric', smooth) - smooth = float(smooth) - if not 0 <= smooth < 1: - raise ValueError('smooth must be between 0 (inclusive) and 1 ' - f'(exclusive), got {smooth}') if image.ndim != 3: raise ValueError(f'3D data must be supplied, got {image.shape}') @@ -1787,8 +1792,7 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None): # force double as passing integer types directly can be problematic! image_shape = image.shape - data_vtk = numpy_support.numpy_to_vtk( - image.ravel(), deep=True, array_type=VTK_DOUBLE) + data_vtk = numpy_to_vtk(image.ravel().astype(float), deep=True) del image mc = vtkDiscreteFlyingEdges3D() @@ -1799,28 +1803,17 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None): imdata.SetOrigin([0, 0, 0]) imdata.GetPointData().SetScalars(data_vtk) - # compute marching cubes + # compute marching cubes on smoothed data mc.SetNumberOfContours(len(level)) for li, lev in enumerate(level): mc.SetValue(li, lev) mc.SetInputData(imdata) - sel_input = mc - if smooth: - smoother = vtkWindowedSincPolyDataFilter() - smoother.SetInputConnection(sel_input.GetOutputPort()) - smoother.SetNumberOfIterations(100) - smoother.BoundarySmoothingOff() - smoother.FeatureEdgeSmoothingOff() - smoother.SetFeatureAngle(120.0) - smoother.SetPassBand(1 - smooth) - smoother.NonManifoldSmoothingOn() - smoother.NormalizeCoordinatesOff() - sel_input = smoother - sel_input.Update() + mc.Update() + mc = _vtk_smooth(mc.GetOutput(), smooth) # get verts and triangles selector = vtkThreshold() - selector.SetInputConnection(sel_input.GetOutputPort()) + selector.SetInputData(mc) dsa = vtkDataSetAttributes() selector.SetInputArrayToProcess( 0, 0, 0, imdata.FIELD_ASSOCIATION_POINTS, dsa.SCALARS) @@ -1839,8 +1832,8 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None): selector.SetUpperThreshold(val) geometry.Update() polydata = geometry.GetOutput() - rr = numpy_support.vtk_to_numpy(polydata.GetPoints().GetData()) - tris = numpy_support.vtk_to_numpy( + rr = vtk_to_numpy(polydata.GetPoints().GetData()) + tris = vtk_to_numpy( polydata.GetPolys().GetConnectivityArray()).reshape(-1, 3) rr = np.ascontiguousarray(rr[:, ::-1]) tris = np.ascontiguousarray(tris[:, ::-1]) @@ -1848,6 +1841,36 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None): return out +def _vtk_smooth(pd, smooth): + _validate_type(smooth, 'numeric', smooth) + smooth = float(smooth) + if not 0 <= smooth < 1: + raise ValueError('smoothing factor must be between 0 (inclusive) and ' + f'1 (exclusive), got {smooth}') + if smooth == 0: + return pd + from vtkmodules.vtkFiltersCore import vtkWindowedSincPolyDataFilter + logger.info(f' Smoothing by a factor of {smooth}') + return_ndarray = False + if isinstance(pd, dict): + pd = _surface_to_polydata(pd) + return_ndarray = True + smoother = vtkWindowedSincPolyDataFilter() + smoother.SetInputData(pd) + smoother.SetNumberOfIterations(100) + smoother.BoundarySmoothingOff() + smoother.FeatureEdgeSmoothingOff() + smoother.SetFeatureAngle(120.0) + smoother.SetPassBand(1 - smooth) + smoother.NonManifoldSmoothingOn() + smoother.NormalizeCoordinatesOff() + smoother.Update() + out = smoother.GetOutput() + if return_ndarray: + out = _polydata_to_surface(out, normals=False) + return out + + def _warn_missing_chs(info, dig_image, after_warp, verbose=None): """Warn that channels are missing.""" # ensure that each electrode contact was marked in at least one voxel diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py index bb60ce17ccc..3a3ffae1cb2 100644 --- a/mne/tests/test_bem.py +++ b/mne/tests/test_bem.py @@ -12,6 +12,7 @@ import pytest from numpy.testing import assert_equal, assert_allclose +import mne from mne import (make_bem_model, read_bem_surfaces, write_bem_surfaces, make_bem_solution, read_bem_solution, write_bem_solution, make_sphere_model, Transform, Info, write_surface, @@ -23,8 +24,9 @@ from mne.utils import catch_logging, check_version from mne.bem import (_ico_downsample, _get_ico_map, _order_surfaces, _assert_complete_surface, _assert_inside, - _check_surface_size, _bem_find_surface) -from mne.surface import read_surface + _check_surface_size, _bem_find_surface, + make_scalp_surfaces) +from mne.surface import read_surface, _get_ico_surface from mne.io import read_info fname_raw = op.join(op.dirname(__file__), '..', 'io', 'tests', 'data', @@ -419,3 +421,52 @@ def test_io_head_bem(tmp_path): assert head['id'] == head_defect['id'] == FIFF.FIFFV_BEM_SURF_ID_HEAD assert np.allclose(head['rr'], head_defect['rr']) assert np.allclose(head['tris'], head_defect['tris']) + + +@pytest.mark.slowtest # ~4 sec locally +def test_make_scalp_surfaces_topology(tmp_path, monkeypatch): + """Test topology checks for make_scalp_surfaces.""" + pytest.importorskip('pyvista') + subjects_dir = tmp_path + subject = 'test' + surf_dir = subjects_dir / subject / 'surf' + makedirs(surf_dir) + surf = _get_ico_surface(2) + surf['rr'] *= 100 # mm + write_surface(surf_dir / 'lh.seghead', surf['rr'], surf['tris']) + + # make it so that decimation really messes up the mesh just by deleting + # the last N tris + def _decimate_surface(points, triangles, n_triangles): + assert len(triangles) >= n_triangles + return points, triangles[:n_triangles] + + monkeypatch.setattr(mne.bem, 'decimate_surface', _decimate_surface) + # TODO: These two errors should probably have the same class... + + # Not enough neighbors + monkeypatch.setattr(mne.bem, '_tri_levels', dict(sparse=315)) + with pytest.raises(ValueError, match='.*have fewer than three.*'): + make_scalp_surfaces(subject, subjects_dir, force=False, verbose=True) + monkeypatch.setattr(mne.bem, '_tri_levels', dict(sparse=319)) + # Incomplete surface (sum of solid angles) + with pytest.raises(RuntimeError, match='.*is not complete.*'): + 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() + + # These are ignorable + monkeypatch.setattr(mne.bem, '_tri_levels', dict(sparse=315)) + with pytest.warns(RuntimeWarning, match='.*have fewer than three.*'): + make_scalp_surfaces( + subject, subjects_dir, force=True, overwrite=True) + surf, = read_bem_surfaces(sparse_path, on_defects='ignore') + assert len(surf['tris']) == 315 + monkeypatch.setattr(mne.bem, '_tri_levels', dict(sparse=319)) + with pytest.warns(RuntimeWarning, match='.*is not complete.*'): + make_scalp_surfaces( + subject, subjects_dir, force=True, overwrite=True) + surf, = read_bem_surfaces(sparse_path, on_defects='ignore') + assert len(surf['tris']) == 319 diff --git a/mne/tests/test_surface.py b/mne/tests/test_surface.py index 7d0913516cc..38ee1f10811 100644 --- a/mne/tests/test_surface.py +++ b/mne/tests/test_surface.py @@ -20,7 +20,7 @@ _normal_orth, _read_patch, _marching_cubes, _voxel_neighbors, warp_montage_volume) from mne.transforms import _get_trans, compute_volume_registration, apply_trans -from mne.utils import (requires_vtk, catch_logging, object_diff, +from mne.utils import (catch_logging, object_diff, requires_freesurfer, requires_nibabel, requires_dipy, _record_warnings) @@ -166,9 +166,9 @@ def test_read_curv(): assert np.logical_or(bin_curv == 0, bin_curv == 1).all() -@requires_vtk def test_decimate_surface_vtk(): """Test triangular surface decimation.""" + pytest.importorskip('pyvista') points = np.array([[-0.00686118, -0.10369860, 0.02615170], [-0.00713948, -0.10370162, 0.02614874], [-0.00686208, -0.10368247, 0.02588313], @@ -225,12 +225,12 @@ def test_normal_orth(): # 0.06 sec locally even with all these params -@requires_vtk @pytest.mark.parametrize('dtype', (np.float64, np.uint16, '>i4')) @pytest.mark.parametrize('value', (1, 12)) @pytest.mark.parametrize('smooth', (0, 0.9)) def test_marching_cubes(dtype, value, smooth): """Test creating surfaces via marching cubes.""" + pytest.importorskip('pyvista') data = np.zeros((50, 50, 50), dtype=dtype) data[20:30, 20:30, 20:30] = value level = [value] diff --git a/mne/utils/__init__.py b/mne/utils/__init__.py index 138eece90eb..4dc09e2cd04 100644 --- a/mne/utils/__init__.py +++ b/mne/utils/__init__.py @@ -45,7 +45,7 @@ requires_good_network, requires_pandas, requires_h5py, ArgvSetter, SilenceStdout, has_freesurfer, has_mne_c, _TempDir, has_nibabel, buggy_mkl_svd, - requires_numpydoc, requires_vtk, requires_freesurfer, + requires_numpydoc, requires_freesurfer, requires_nitime, requires_dipy, requires_neuromag2ft, requires_pylsl, assert_object_equal, assert_and_remove_boundary_annot, diff --git a/mne/utils/_testing.py b/mne/utils/_testing.py index 33ace6818ca..b3068d3d8b3 100644 --- a/mne/utils/_testing.py +++ b/mne/utils/_testing.py @@ -91,7 +91,7 @@ def requires_version(library, min_version='0.0'): import pytest reason = f'Requires {library}' if min_version != '0.0': - reason += ' version >= {min_version}' + reason += f' version >= {min_version}' return pytest.mark.skipif(not check_version(library, min_version), reason=reason) @@ -153,7 +153,6 @@ def requires_freesurfer(arg): requires_neuromag2ft = partial(requires_module, name='neuromag2ft', call=_n2ft_call) -requires_vtk = partial(requires_module, name='vtk') requires_good_network = partial( requires_module, name='good network connection', call='if int(os.environ.get("MNE_SKIP_NETWORK_TESTS", 0)):\n' diff --git a/mne/utils/misc.py b/mne/utils/misc.py index 5a071b1b89c..af137651f20 100644 --- a/mne/utils/misc.py +++ b/mne/utils/misc.py @@ -24,10 +24,10 @@ from ._logging import logger, verbose, warn -def _pl(x, non_pl=''): +def _pl(x, non_pl='', pl='s'): """Determine if plural should be used.""" len_x = x if isinstance(x, (int, np.generic)) else len(x) - return non_pl if len_x == 1 else 's' + return non_pl if len_x == 1 else pl def _explain_exception(start=-1, stop=None, prefix='> '): diff --git a/mne/viz/backends/_pyvista.py b/mne/viz/backends/_pyvista.py index 311b2b17bda..3c2c81bafd1 100644 --- a/mne/viz/backends/_pyvista.py +++ b/mne/viz/backends/_pyvista.py @@ -18,7 +18,6 @@ import warnings import numpy as np -import vtk from ._abstract import _AbstractRenderer, Figure3D from ._utils import (_get_colormap_from_array, _alpha_blend_background, @@ -38,7 +37,30 @@ except ImportError: from pyvista import BackgroundPlotter from pyvista.plotting.plotting import _ALL_PLOTTERS -VTK9 = _compare_version(getattr(vtk, 'VTK_VERSION', '9.0'), '>=', '9.0') + +from vtkmodules.vtkCommonCore import ( + vtkCommand, vtkLookupTable, VTK_UNSIGNED_CHAR) +from vtkmodules.vtkCommonDataModel import VTK_VERTEX, vtkPiecewiseFunction +from vtkmodules.vtkCommonTransforms import vtkTransform +from vtkmodules.vtkFiltersCore import vtkCellDataToPointData, vtkGlyph3D +from vtkmodules.vtkFiltersGeneral import ( + vtkTransformPolyDataFilter, vtkMarchingContourFilter) +from vtkmodules.vtkFiltersHybrid import vtkPolyDataSilhouette +from vtkmodules.vtkFiltersSources import ( + vtkSphereSource, vtkConeSource, vtkCylinderSource, vtkArrowSource, + vtkPlatonicSolidSource, vtkGlyphSource2D) +from vtkmodules.vtkImagingCore import vtkImageReslice +from vtkmodules.vtkRenderingCore import ( + vtkMapper, vtkActor, vtkCellPicker, vtkColorTransferFunction, + vtkPolyDataMapper, vtkVolume, vtkCoordinate, vtkDataSetMapper) +from vtkmodules.vtkRenderingVolumeOpenGL2 import vtkSmartVolumeMapper +from vtkmodules.util.numpy_support import numpy_to_vtk +try: + from vtkmodules.vtkCommonCore import VTK_VERSION +except Exception: # some bad versions of VTK + VTK_VERSION = '9.0' + +VTK9 = _compare_version(VTK_VERSION, '>=', '9.0') _FIGURES = dict() @@ -278,7 +300,9 @@ def set_interaction(self, interaction): if hasattr(self.plotter, 'enable_rubber_band_2d_style'): self.plotter.enable_rubber_band_2d_style() else: - style = vtk.vtkInteractorStyleRubberBand2D() + from vtkmodules.vtkInteractionStyle import\ + vtkInteractorStyleRubberBand2D + style = vtkInteractorStyleRubberBand2D() self.plotter.interactor.SetInteractorStyle(style) else: for renderer in self._all_renderers: @@ -427,6 +451,7 @@ def surface(self, surface, color=None, opacity=1.0, def sphere(self, center, color, scale, opacity=1.0, resolution=8, backface_culling=False, radius=None): + from vtkmodules.vtkFiltersSources import vtkSphereSource factor = 1.0 if radius is not None else scale center = np.array(center, dtype=float) if len(center) == 0: @@ -435,7 +460,7 @@ def sphere(self, center, color, scale, opacity=1.0, _check_option('center.shape[-1]', center.shape[-1], (3,)) with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=FutureWarning) - sphere = vtk.vtkSphereSource() + sphere = vtkSphereSource() sphere.SetThetaResolution(resolution) sphere.SetPhiResolution(resolution) if radius is not None: @@ -494,7 +519,7 @@ def quiver3d(self, x, y, z, u, v, w, color, scale, mode, resolution=8, vectors = np.c_[u, v, w] points = np.vstack(np.c_[x, y, z]) n_points = len(points) - cell_type = np.full(n_points, vtk.VTK_VERTEX) + cell_type = np.full(n_points, VTK_VERTEX) cells = np.c_[np.full(n_points, 1), range(n_points)] args = (cells, cell_type, points) if not VTK9: @@ -517,20 +542,20 @@ def quiver3d(self, x, y, z, u, v, w, color, scale, mode, resolution=8, else: tr = None if mode == 'cone': - glyph = vtk.vtkConeSource() + glyph = vtkConeSource() glyph.SetCenter(0.5, 0, 0) if glyph_radius is not None: glyph.SetRadius(glyph_radius) elif mode == 'cylinder': - glyph = vtk.vtkCylinderSource() + glyph = vtkCylinderSource() if glyph_radius is not None: glyph.SetRadius(glyph_radius) elif mode == 'oct': - glyph = vtk.vtkPlatonicSolidSource() + glyph = vtkPlatonicSolidSource() glyph.SetSolidTypeToOctahedron() else: assert mode == 'sphere', mode # guaranteed above - glyph = vtk.vtkSphereSource() + glyph = vtkSphereSource() if mode == 'cylinder': if glyph_height is not None: glyph.SetHeight(glyph_height) @@ -538,18 +563,18 @@ def quiver3d(self, x, y, z, u, v, w, color, scale, mode, resolution=8, glyph.SetCenter(glyph_center) if glyph_resolution is not None: glyph.SetResolution(glyph_resolution) - tr = vtk.vtkTransform() + tr = vtkTransform() tr.RotateWXYZ(90, 0, 0, 1) elif mode == 'oct': if solid_transform is not None: assert solid_transform.shape == (4, 4) - tr = vtk.vtkTransform() + tr = vtkTransform() tr.SetMatrix( solid_transform.astype(np.float64).ravel()) if tr is not None: # fix orientation glyph.Update() - trp = vtk.vtkTransformPolyDataFilter() + trp = vtkTransformPolyDataFilter() trp.SetInputData(glyph.GetOutput()) trp.SetTransform(tr) glyph = trp @@ -614,9 +639,9 @@ def text3d(self, x, y, z, text, scale, color='white'): def scalarbar(self, source, color="white", title=None, n_labels=4, bgcolor=None, **extra_kwargs): - if isinstance(source, vtk.vtkMapper): + if isinstance(source, vtkMapper): mapper = source - elif isinstance(source, vtk.vtkActor): + elif isinstance(source, vtkActor): mapper = source.GetMapper() else: mapper = None @@ -705,7 +730,7 @@ def _disabled_interaction(self): self.plotter.enable() def _actor(self, mapper=None): - actor = vtk.vtkActor() + actor = vtkActor() if mapper is not None: actor.SetMapper(mapper) _hide_testing_actor(actor) @@ -721,12 +746,12 @@ def _update_picking_callback(self, on_button_release, on_pick): add_obs = self.plotter.iren.add_observer - add_obs(vtk.vtkCommand.RenderEvent, on_mouse_move) - add_obs(vtk.vtkCommand.LeftButtonPressEvent, on_button_press) - add_obs(vtk.vtkCommand.EndInteractionEvent, on_button_release) - self.plotter.picker = vtk.vtkCellPicker() + add_obs(vtkCommand.RenderEvent, on_mouse_move) + add_obs(vtkCommand.LeftButtonPressEvent, on_button_press) + add_obs(vtkCommand.EndInteractionEvent, on_button_release) + self.plotter.picker = vtkCellPicker() self.plotter.picker.AddObserver( - vtk.vtkCommand.EndPickEvent, + vtkCommand.EndPickEvent, on_pick ) self.plotter.picker.SetVolumeOpacityIsovalue(0.) @@ -740,7 +765,6 @@ def _set_mesh_scalars(self, mesh, scalars, name): def _set_colormap_range(self, actor, ctable, scalar_bar, rng=None, background_color=None): - from vtk.util.numpy_support import numpy_to_vtk if rng is not None: mapper = actor.GetMapper() mapper.SetScalarRange(*rng) @@ -751,15 +775,12 @@ def _set_colormap_range(self, actor, ctable, scalar_bar, rng=None, if background_color is not None: background_color = np.array(background_color) * 255 ctable = _alpha_blend_background(ctable, background_color) - lut.SetTable(numpy_to_vtk(ctable, - array_type=vtk.VTK_UNSIGNED_CHAR)) + lut.SetTable(numpy_to_vtk(ctable, array_type=VTK_UNSIGNED_CHAR)) lut.SetRange(*rng) def _set_volume_range(self, volume, ctable, alpha, scalar_bar, rng): - import vtk - from vtk.util.numpy_support import numpy_to_vtk - color_tf = vtk.vtkColorTransferFunction() - opacity_tf = vtk.vtkPiecewiseFunction() + color_tf = vtkColorTransferFunction() + opacity_tf = vtkPiecewiseFunction() for loc, color in zip(np.linspace(*rng, num=len(ctable)), ctable): color_tf.AddRGBPoint(loc, *(color[:-1] / 255.)) opacity_tf.AddPoint(loc, color[-1] * alpha / 255.) @@ -771,13 +792,14 @@ def _set_volume_range(self, volume, ctable, alpha, scalar_bar, rng): prop.ShadeOn() prop.SetInterpolationTypeToLinear() if scalar_bar is not None: - lut = vtk.vtkLookupTable() + lut = vtkLookupTable() lut.SetRange(*rng) lut.SetTable(numpy_to_vtk(ctable)) scalar_bar.SetLookupTable(lut) def _sphere(self, center, color, radius): - sphere = vtk.vtkSphereSource() + from vtkmodules.vtkFiltersSources import vtkSphereSource + sphere = vtkSphereSource() sphere.SetThetaResolution(8) sphere.SetPhiResolution(8) sphere.SetRadius(radius) @@ -803,7 +825,7 @@ def _volume(self, dimensions, origin, spacing, scalars, # Add contour of enclosed volume (use GetOutput instead of # GetOutputPort below to avoid updating) if surface_alpha > 0 or resolution is not None: - grid_alg = vtk.vtkCellDataToPointData() + grid_alg = vtkCellDataToPointData() grid_alg.SetInputDataObject(grid) grid_alg.SetPassCellData(False) grid_alg.Update() @@ -811,23 +833,23 @@ def _volume(self, dimensions, origin, spacing, scalars, grid_alg = None if surface_alpha > 0: - grid_surface = vtk.vtkMarchingContourFilter() + grid_surface = vtkMarchingContourFilter() grid_surface.ComputeNormalsOn() grid_surface.ComputeScalarsOff() grid_surface.SetInputData(grid_alg.GetOutput()) grid_surface.SetValue(0, 0.1) grid_surface.Update() - grid_mesh = vtk.vtkPolyDataMapper() + grid_mesh = vtkPolyDataMapper() grid_mesh.SetInputData(grid_surface.GetOutput()) else: grid_mesh = None - mapper = vtk.vtkSmartVolumeMapper() + mapper = vtkSmartVolumeMapper() if resolution is None: # native mapper.SetScalarModeToUseCellData() mapper.SetInputDataObject(grid) else: - upsampler = vtk.vtkImageReslice() + upsampler = vtkImageReslice() upsampler.SetInterpolationModeToLinear() # default anyway upsampler.SetOutputSpacing(*([resolution] * 3)) upsampler.SetInputConnection(grid_alg.GetOutputPort()) @@ -835,20 +857,20 @@ def _volume(self, dimensions, origin, spacing, scalars, # Additive, AverageIntensity, and Composite might also be reasonable remap = dict(composite='Composite', mip='MaximumIntensity') getattr(mapper, f'SetBlendModeTo{remap[blending]}')() - volume_pos = vtk.vtkVolume() + volume_pos = vtkVolume() volume_pos.SetMapper(mapper) dist = grid.length / (np.mean(grid.dimensions) - 1) volume_pos.GetProperty().SetScalarOpacityUnitDistance(dist) if center is not None and blending == 'mip': # We need to create a minimum intensity projection for the neg half - mapper_neg = vtk.vtkSmartVolumeMapper() + mapper_neg = vtkSmartVolumeMapper() if resolution is None: # native mapper_neg.SetScalarModeToUseCellData() mapper_neg.SetInputDataObject(grid) else: mapper_neg.SetInputConnection(upsampler.GetOutputPort()) mapper_neg.SetBlendModeToMinimumIntensity() - volume_neg = vtk.vtkVolume() + volume_neg = vtkVolume() volume_neg.SetMapper(mapper_neg) volume_neg.GetProperty().SetScalarOpacityUnitDistance(dist) else: @@ -858,11 +880,11 @@ def _volume(self, dimensions, origin, spacing, scalars, def _silhouette(self, mesh, color=None, line_width=None, alpha=None, decimate=None): mesh = mesh.decimate(decimate) if decimate is not None else mesh - silhouette_filter = vtk.vtkPolyDataSilhouette() + silhouette_filter = vtkPolyDataSilhouette() silhouette_filter.SetInputData(mesh) silhouette_filter.SetCamera(self.plotter.renderer.GetActiveCamera()) silhouette_filter.SetEnableFeatureAngle(0) - silhouette_mapper = vtk.vtkPolyDataMapper() + silhouette_mapper = vtkPolyDataMapper() silhouette_mapper.SetInputConnection( silhouette_filter.GetOutputPort()) actor, prop = self.plotter.add_actor( @@ -940,8 +962,7 @@ def _mat_to_array(vtk_mat): def _3d_to_2d(plotter, xyz): # https://vtk.org/Wiki/VTK/Examples/Cxx/Utilities/Coordinate - import vtk - coordinate = vtk.vtkCoordinate() + coordinate = vtkCoordinate() coordinate.SetCoordinateSystemToWorld() xy = list() for coord in xyz: @@ -1078,19 +1099,19 @@ def _process_events(plotter): def _add_camera_callback(camera, callback): - camera.AddObserver(vtk.vtkCommand.ModifiedEvent, callback) + camera.AddObserver(vtkCommand.ModifiedEvent, callback) def _arrow_glyph(grid, factor): - glyph = vtk.vtkGlyphSource2D() + glyph = vtkGlyphSource2D() glyph.SetGlyphTypeToArrow() glyph.FilledOff() glyph.Update() # fix position - tr = vtk.vtkTransform() + tr = vtkTransform() tr.Translate(0.5, 0., 0.) - trp = vtk.vtkTransformPolyDataFilter() + trp = vtkTransformPolyDataFilter() trp.SetInputConnection(glyph.GetOutputPort()) trp.SetTransform(tr) trp.Update() @@ -1103,7 +1124,7 @@ def _arrow_glyph(grid, factor): factor=factor, geom=trp.GetOutputPort(), ) - mapper = vtk.vtkDataSetMapper() + mapper = vtkDataSetMapper() mapper.SetInputConnection(alg.GetOutputPort()) return mapper @@ -1111,10 +1132,10 @@ def _arrow_glyph(grid, factor): def _glyph(dataset, scale_mode='scalar', orient=True, scalars=True, factor=1.0, geom=None, tolerance=0.0, absolute=False, clamping=False, rng=None): if geom is None: - arrow = vtk.vtkArrowSource() + arrow = vtkArrowSource() arrow.Update() geom = arrow.GetOutputPort() - alg = vtk.vtkGlyph3D() + alg = vtkGlyph3D() alg.SetSourceConnection(geom) if isinstance(scalars, str): dataset.active_scalars_name = scalars diff --git a/tools/azure_dependencies.sh b/tools/azure_dependencies.sh index 0f5c855c40a..34b2b2ff0c9 100755 --- a/tools/azure_dependencies.sh +++ b/tools/azure_dependencies.sh @@ -20,6 +20,7 @@ elif [ "${TEST_MODE}" == "pip-pre" ]; then python -m pip install --progress-bar off --upgrade --pre --only-binary ":all:" https://github.com/pyvista/pyvista-wheels/raw/3b70f3933bc246b035354b172a0459ffc96b0343/vtk-9.1.0.dev0-cp310-cp310-win_amd64.whl python -m pip install --progress-bar off https://github.com/pyvista/pyvista/zipball/main python -m pip install --progress-bar off https://github.com/pyvista/pyvistaqt/zipball/main + python -m pip install --progress-bar off imageio-ffmpeg xlrd mffpy python-picard patsy EXTRA_ARGS="--pre" else echo "Unknown run type ${TEST_MODE}" diff --git a/tools/circleci_dependencies.sh b/tools/circleci_dependencies.sh index 564f0576ad0..b6d6c61b5c1 100755 --- a/tools/circleci_dependencies.sh +++ b/tools/circleci_dependencies.sh @@ -15,6 +15,6 @@ else # standard doc build # TODO: Revert to upstream/main once https://github.com/mne-tools/mne-qt-browser/pull/105 is merged python -m pip install --upgrade --progress-bar off https://github.com/mne-tools/mne-qt-browser/zipball/main https://github.com/sphinx-gallery/sphinx-gallery/zipball/master https://github.com/pyvista/pyvista/zipball/main https://github.com/pyvista/pyvistaqt/zipball/main # deal with comparisons and escapes (https://app.circleci.com/pipelines/github/mne-tools/mne-python/9686/workflows/3fd32b47-3254-4812-8b9a-8bab0d646d18/jobs/32934) - python -m pip install --upgrade quantities + python -m pip install --upgrade --progress-bar off quantities fi python -m pip install -e . diff --git a/tools/github_actions_dependencies.sh b/tools/github_actions_dependencies.sh index 39e6850176c..f56f7119a52 100755 --- a/tools/github_actions_dependencies.sh +++ b/tools/github_actions_dependencies.sh @@ -32,8 +32,8 @@ else pip install --progress-bar off https://github.com/pyvista/pyvista/zipball/main echo "pyvistaqt" pip install --progress-bar off https://github.com/pyvista/pyvistaqt/zipball/main - echo "imageio-ffmpeg, xlrd, mffpy" - pip install --progress-bar off --pre imageio-ffmpeg xlrd mffpy + echo "imageio-ffmpeg, xlrd, mffpy, python-picard" + pip install --progress-bar off --pre imageio-ffmpeg xlrd mffpy python-picard patsy if [ "$OSTYPE" == "darwin"* ]; then echo "pyobjc-framework-Cocoa" pip install --progress-bar off pyobjc-framework-Cocoa>=5.2.0 diff --git a/tutorials/forward/80_fix_bem_in_blender.py b/tutorials/forward/80_fix_bem_in_blender.py index d00a33b9ece..b319ba74bd6 100644 --- a/tutorials/forward/80_fix_bem_in_blender.py +++ b/tutorials/forward/80_fix_bem_in_blender.py @@ -2,22 +2,23 @@ """ .. _tut-fix-meshes: -=============================== -Editing BEM surfaces in Blender -=============================== +============================ +Fixing BEM and head surfaces +============================ Sometimes when creating a BEM model the surfaces need manual correction because of a series of problems that can arise (e.g. intersection between surfaces). Here, we will see how this can be achieved by exporting the surfaces to the 3D modeling program `Blender `_, editing them, and -re-importing them. +re-importing them. We will also give a simple example of how to use +:ref:`pymeshfix ` to fix topological problems. -This tutorial is based on https://github.com/ezemikulan/blender_freesurfer by -Ezequiel Mikulan. +Much of this tutorial is based on +https://github.com/ezemikulan/blender_freesurfer by Ezequiel Mikulan. .. contents:: Page contents :local: - :depth: 2 + :depth: 3 """ @@ -32,13 +33,14 @@ # sphinx_gallery_thumbnail_path = '_static/blender_import_obj/blender_import_obj2.jpg' # noqa import os -import os.path as op import shutil import mne data_path = mne.datasets.sample.data_path() -subjects_dir = op.join(data_path, 'subjects') -bem_dir = op.join(subjects_dir, 'sample', 'bem', 'flash') +subjects_dir = data_path / 'subjects' +bem_dir = subjects_dir / 'sample' / 'bem' / 'flash' +surf_dir = subjects_dir / 'sample' / 'surf' + # %% # Exporting surfaces to Blender # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -53,24 +55,22 @@ # folder called ``conv`` inside the FreeSurfer subject folder to keep them in. # Put the converted surfaces in a separate 'conv' folder -conv_dir = op.join(subjects_dir, 'sample', 'conv') +conv_dir = subjects_dir / 'sample' / 'conv' os.makedirs(conv_dir, exist_ok=True) # Load the inner skull surface and create a problem # The metadata is empty in this example. In real study, we want to write the # original metadata to the fixed surface file. Set read_metadata=True to do so. -coords, faces = mne.read_surface(op.join(bem_dir, 'inner_skull.surf')) +coords, faces = mne.read_surface(bem_dir / 'inner_skull.surf') coords[0] *= 1.1 # Move the first vertex outside the skull # Write the inner skull surface as an .obj file that can be imported by # Blender. -mne.write_surface(op.join(conv_dir, 'inner_skull.obj'), coords, faces, - overwrite=True) +mne.write_surface(conv_dir / 'inner_skull.obj', coords, faces, overwrite=True) # Also convert the outer skull surface. -coords, faces = mne.read_surface(op.join(bem_dir, 'outer_skull.surf')) -mne.write_surface(op.join(conv_dir, 'outer_skull.obj'), coords, faces, - overwrite=True) +coords, faces = mne.read_surface(bem_dir / 'outer_skull.surf') +mne.write_surface(conv_dir / 'outer_skull.obj', coords, faces, overwrite=True) # %% # Editing in Blender @@ -122,9 +122,9 @@ # In order to be able to run this tutorial script top to bottom, we here # simulate the edits you did manually in Blender using Python code: -coords, faces = mne.read_surface(op.join(conv_dir, 'inner_skull.obj')) +coords, faces = mne.read_surface(conv_dir / 'inner_skull.obj') coords[0] /= 1.1 # Move the first vertex back inside the skull -mne.write_surface(op.join(conv_dir, 'inner_skull_fixed.obj'), coords, faces, +mne.write_surface(conv_dir / 'inner_skull_fixed.obj', coords, faces, overwrite=True) # %% @@ -135,19 +135,18 @@ # surfaces in case you make a mistake! # Read the fixed surface -coords, faces = mne.read_surface(op.join(conv_dir, 'inner_skull_fixed.obj')) +coords, faces = mne.read_surface(conv_dir / 'inner_skull_fixed.obj') # Backup the original surface -shutil.copy(op.join(bem_dir, 'inner_skull.surf'), - op.join(bem_dir, 'inner_skull_orig.surf')) +shutil.copy(bem_dir / 'inner_skull.surf', bem_dir / 'inner_skull_orig.surf') # Overwrite the original surface with the fixed version # In real study you should provide the correct metadata using ``volume_info=`` # This could be accomplished for example with: # -# _, _, vol_info = mne.read_surface(op.join(bem_dir, 'inner_skull.surf'), +# _, _, vol_info = mne.read_surface(bem_dir / 'inner_skull.surf', # read_metadata=True) -# mne.write_surface(op.join(bem_dir, 'inner_skull.surf'), coords, faces, +# mne.write_surface(bem_dir / 'inner_skull.surf', coords, faces, # volume_info=vol_info, overwrite=True) # %% @@ -165,15 +164,15 @@ # ``-head.fif`` from the edited surface file for coregistration. # Load the fixed surface -coords, faces = mne.read_surface(op.join(bem_dir, 'outer_skin.surf')) +coords, faces = mne.read_surface(bem_dir / 'outer_skin.surf') # Make sure we are in the correct directory -head_dir = op.dirname(bem_dir) +head_dir = bem_dir.parent # Remember to backup the original head file in advance! # Overwrite the original head file # -# mne.write_head_bem(op.join(head_dir, 'sample-head.fif'), coords, faces, +# mne.write_head_bem(head_dir / 'sample-head.fif', coords, faces, # overwrite=True) # %% @@ -187,29 +186,33 @@ # If ``-head-dense.fif`` does not exist, you need to run # ``mne make_scalp_surfaces`` first. # [0] because a list of surfaces is returned -surf = mne.read_bem_surfaces(op.join(head_dir, 'sample-head.fif'))[0] +surf = mne.read_bem_surfaces(head_dir / 'sample-head.fif')[0] # For consistency only coords = surf['rr'] faces = surf['tris'] # Write the head as an .obj file for editing -mne.write_surface(op.join(conv_dir, 'sample-head.obj'), +mne.write_surface(conv_dir / 'sample-head.obj', coords, faces, overwrite=True) # Usually here you would go and edit your meshes. # # Here we just use the same surface as if it were fixed # Read in the .obj file -coords, faces = mne.read_surface(op.join(conv_dir, 'sample-head.obj')) +coords, faces = mne.read_surface(conv_dir / 'sample-head.obj') # Remember to backup the original head file in advance! # Overwrite the original head file # -# mne.write_head_bem(op.join(head_dir, 'sample-head.fif'), coords, faces, +# mne.write_head_bem(head_dir / 'sample-head.fif', coords, faces, # overwrite=True) # %% +# .. note:: See also :ref:`tut-fix-meshes-smoothing` for a possible alternative +# high-resolution head fix using FreeSurfer smoothing instead of +# blender. +# # Blender editing tips # ~~~~~~~~~~~~~~~~~~~~ # @@ -217,85 +220,136 @@ # restrict one surface inside another surface (for example the Skull inside the # Outer Skin). Here is how to use it: # -# (1) Select the surface that is creating the problem. (2) In *Edit Mode*, -# press :kbd:`C` to use the circle selection tool to select the vertices that -# are outside. (3-5) In the *Object Data Properties* tab use the ``+`` button -# to add a *Vertex Group* and click *Assign* to assign the current selection to -# the group. (6-8) In the *Modifiers* tab go to *Add Modifier* add a -# *Shrinkwrap* modifier and set it to snap *Inside* with the outer surface as -# the *Target* and the *Group* that you created before as the *Vertex Group*. -# You can then use the *Offset* parameter to adjust the distance. (9) In -# *Object Mode* click on the down-pointing arrow of the *Shrinkwrap* modifier -# and click on *Apply*. +# ① Select surface +# Select the surface that is creating the problem. +# ② Select vertices +# In *Edit Mode*, press :kbd:`C` to use the circle selection tool to +# select the vertices that are outside. +# ③-⑤ Group vertices +# In the *Object Data Properties* tab use the ``+`` button +# to add a *Vertex Group* and click *Assign* to assign the current +# selection to the group. +# ⑥-⑧ Shrinkwrap +# In the *Modifiers* tab go to *Add Modifier* add a +# *Shrinkwrap* modifier and set it to snap *Inside* with the outer surface +# as the *Target* and the *Group* that you created before as the +# *Vertex Group*. You can then use the *Offset* parameter to adjust the +# distance. +# ⑨ Apply +# In *Object Mode* click on the down-pointing arrow of the *Shrinkwrap* +# modifier and click on *Apply*. # # .. image:: ../../_static/blender_import_obj/blender_import_obj4.jpg # :alt: Shrinkwrap functionality in Blender - -# %% +# # That's it! You are ready to continue with your analysis pipeline (e.g. # running :func:`mne.make_bem_model`). # # What if you still get an error? -# --------------------------------- +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # +# Blender operations +# ~~~~~~~~~~~~~~~~~~ # When editing BEM surfaces/meshes in Blender, make sure to use # tools that do not change the number or order of vertices, or the geometry # of triangular faces. For example, avoid the extrusion tool, because it # duplicates the extruded vertices. # -# Below are some examples of errors you might encounter when running the -# `mne.make_bem_model` function, and the likely causes of those errors. +# .. _tut-fix-meshes-smoothing: # +# Cleaning up a bad dense head surface by smoothing +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# If you get a rough head surfaces when using :ref:`mne make_scalp_surfaces`, +# consider smoothing your T1 ahead of time with a Gaussian kernel with +# FreeSurfer using something like the following within the subject's ``mri`` +# directory: # -# 1. Cannot decimate to requested ico grade +# .. code-block:: console # -# This error is caused by having too few or too many vertices. The full -# error is something like: +# $ mri_convert --fwhm 3 T1.mgz T1_smoothed_3.mgz # -# .. code-block:: console +# Here the ``--fwhm`` argument determines how much smoothing (in mm) to apply. +# Then delete ``SUBJECTS_DIR/SUBJECT/surf/lh.seghead``, and re-run +# :ref:`mne make_scalp_surfaces` with the additional arguments +# ``--mri="T1_smoothed_3.mgz" --overwrite``, and you should get cleaner +# surfaces. # -# RuntimeError: Cannot decimate to requested ico grade 4. The provided -# BEM surface has 20516 triangles, which cannot be isomorphic with a -# subdivided icosahedron. Consider manually decimating the surface to a -# suitable density and then use ico=None in make_bem_model. +# Topological errors +# ~~~~~~~~~~~~~~~~~~ +# MNE-Python requires that meshes satisfy some topological checks to ensure +# that subsequent processing like BEM solving and electrode projection can +# work properly. # -# 2. Surface inner skull has topological defects +# Below are some examples of errors you might encounter when working with +# meshes in MNE-Python, and the likely causes of those errors. # -# This error can occur when trying to match the original number of -# triangles by removing vertices. The full error looks like: +# 1. Cannot decimate to requested ico grade +# This error is caused by having too few or too many vertices. The full +# error is something like: # -# .. code-block:: console +# .. code-block:: text # -# RuntimeError: Surface inner skull has topological defects: 12 / 20484 -# vertices have fewer than three neighboring triangles [733, 1014, 2068, -# 7732, 8435, 8489, 10181, 11120, 11121, 11122, 11304, 11788] +# RuntimeError: Cannot decimate to requested ico grade 4. The provided +# BEM surface has 20516 triangles, which cannot be isomorphic with a +# subdivided icosahedron. Consider manually decimating the surface to a +# suitable density and then use ico=None in make_bem_model. # -# 3. Surface inner skull is not complete +# 2. Surface inner skull has topological defects +# This error can occur when trying to match the original number of +# triangles by removing vertices. The full error looks like: # -# This error (like the previous error) reflects a problem with the surface -# topology (i.e., the expected pattern of vertices/edges/faces is -# disrupted). +# .. code-block:: text # -# .. code-block:: console +# RuntimeError: Surface inner skull has topological defects: 12 / 20484 +# vertices have fewer than three neighboring triangles [733, 1014, +# 2068, 7732, 8435, 8489, 10181, 11120, 11121, 11122, 11304, 11788] # -# RuntimeError: Surface inner skull is not complete (sum of solid -# angles yielded 0.999668, should be 1.) +# 3. Surface inner skull is not complete +# This error (like the previous error) reflects a problem with the surface +# topology (i.e., the expected pattern of vertices/edges/faces is +# disrupted). # -# 4. Triangle ordering is wrong +# .. code-block:: text # -# This error reflects a mismatch between how the surface is represented in -# memory (the order of the vertex/face definitions) and what is expected by -# MNE-Python. The full error is: +# RuntimeError: Surface inner skull is not complete (sum of solid +# angles yielded 0.999668, should be 1.) # -# .. code-block:: console +# 4. Triangle ordering is wrong +# This error reflects a mismatch between how the surface is represented in +# memory (the order of the vertex/face definitions) and what is expected +# by MNE-Python. The full error is: # -# RuntimeError: The source surface has a matching number of -# triangles but ordering is wrong +# .. code-block:: text # +# RuntimeError: The source surface has a matching number of triangles +# but ordering is wrong # +# Dealing with topology in blender +# -------------------------------- # For any of these errors, it is usually easiest to start over with the # unedited BEM surface and try again, making sure to only *move* vertices and # faces without *adding* or *deleting* any. For example, # select a circle of vertices, then press :kbd:`G` to drag them to the desired # location. Smoothing a group of selected vertices in Blender (by # right-clicking and selecting "Smooth Vertices") can also be helpful. +# +# .. _tut-fix-meshes-pymeshfix: +# +# Dealing with topology using pymeshfix +# ------------------------------------- +# `pymeshfix `__ is a GPL-licensed Python +# module designed to produce water-tight meshes that satisfy the topological +# checks listed above. For example, if your +# ``'SUBJECTS_DIR/SUBJECT/surf/lh.seghead'`` has topological problems and +# :ref:`tut-fix-meshes-smoothing` does not work, you can try fixing the mesh +# instead. After installing ``pymeshfix`` using conda or pip, from within the +# ``'SUBJECTS_DIR/SUBJECT/surf'`` directory, you could try:: +# +# >>> import pymeshfix +# >>> rr, tris = mne.read_surface('lh.seghead') +# >>> rr, tris = pymeshfix.clean_from_arrays(rr, tris) +# >>> mne.write_surface('lh.seghead', rr, tris) +# +# In some cases this could fix the topology such that a subsequent call to +# :ref:`mne make_scalp_surfaces` will succeed without needing to use the +# ``--force`` parameter.