diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 78580fe0ae0..9c4ddef1613 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -53,6 +53,8 @@ Changelog - Add function to check the type of a FIF file using :func:`mne.what` and `mne what ` by `Eric Larson`_ +- Add support for specifying the initial time and/or position and providing a :class:`mne.SourceMorph` instead of :class:`mne.SourceSpaces` in :func:`mne.viz.plot_volume_source_estimates` by `Eric Larson`_ + - Speed up morph map generation in :func:`mne.read_morph_map` by ~5-10x by using :func:`numba.jit` by `Eric Larson`_ - Speed up :func:`mne.setup_volume_source_space`, especially when ``volume_label is not None`` by `Eric Larson`_ @@ -90,6 +92,14 @@ Bug - Fix one-sample baseline issue in :class:`mne.BaseEpochs` when using `tmin=0` by `Milan Rybář`_ +- Fix bug in :func:`mne.viz.plot_volume_source_estimates` where ``'glass_brain'`` MRIs were not transformed to MNI space, by `Eric Larson`_ + +- Fix bug in :func:`mne.viz.plot_volume_source_estimates` where MRIs with voxels not in RAS orientation could not be browsed properly, by `Eric Larson`_ + +- Fix bug in :meth:`mne.SourceMorph.apply` where output STCs had ``stc.vertices`` defined improperly, by `Eric Larson`_ + +- Fix bug in :meth:`mne.SourceMorph.apply` where the default was errantly ``mri_space=False`` instead of ``mri_space=None`` (as documented), by `Eric Larson`_ + - Fix :meth:`mne.io.Raw.set_annotations` for ``meas_date`` previous to 1970 by `Joan Massich`_ - Fix horizontal spacing issues in :meth:`mne.io.Raw.plot_psd` by `Jeff Hanna`_ diff --git a/doc/conf.py b/doc/conf.py index a6dfadf63fa..f90dd631ced 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -533,9 +533,11 @@ def reset_warnings(gallery_conf, fname): # nibabel 'Nifti1Image': 'nibabel.nifti1.Nifti1Image', 'Nifti2Image': 'nibabel.nifti2.Nifti2Image', + 'SpatialImage': 'nibabel.spatialimages.SpatialImage', # MNE 'Label': 'mne.Label', 'Forward': 'mne.Forward', 'Evoked': 'mne.Evoked', 'Info': 'mne.Info', 'SourceSpaces': 'mne.SourceSpaces', + 'SourceMorph': 'mne.SourceMorph', 'Epochs': 'mne.Epochs', 'Layout': 'mne.channels.Layout', 'EvokedArray': 'mne.EvokedArray', 'BiHemiLabel': 'mne.BiHemiLabel', 'AverageTFR': 'mne.time_frequency.AverageTFR', diff --git a/examples/inverse/plot_lcmv_beamformer_volume.py b/examples/inverse/plot_lcmv_beamformer_volume.py index fe873e7180e..612118bd6d4 100644 --- a/examples/inverse/plot_lcmv_beamformer_volume.py +++ b/examples/inverse/plot_lcmv_beamformer_volume.py @@ -3,14 +3,13 @@ Compute LCMV inverse solution in volume source space ==================================================== -Compute LCMV beamformer on an auditory evoked dataset in a volume source space. +Compute LCMV beamformer on an auditory evoked dataset in a volume source space, +and show activation on ``fsaverage``. """ # Author: Alexandre Gramfort # # License: BSD (3-clause) -# sphinx_gallery_thumbnail_number = 3 - import mne from mne.datasets import sample from mne.beamformer import make_lcmv, apply_lcmv @@ -22,8 +21,7 @@ data_path = sample.data_path() subjects_dir = data_path + '/subjects' -raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' -event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif' +raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-vol-7-fwd.fif' # Get epochs @@ -35,7 +33,7 @@ # Setup for reading the raw data raw = mne.io.read_raw_fif(raw_fname, preload=True) raw.info['bads'] = ['MEG 2443', 'EEG 053'] # 2 bads channels -events = mne.read_events(event_fname) +events = mne.find_events(raw) # Pick the channels of interest raw.pick(['meg', 'eog']) @@ -83,15 +81,34 @@ # You can save result in stc files with: # stc.save('lcmv-vol') - -clim = dict(kind='value', pos_lims=[0.3, 0.6, 0.9]) -stc.plot(src=forward['src'], subject='sample', subjects_dir=subjects_dir, - clim=clim) +lims = [0.3, 0.6, 0.9] +stc.plot( + src=forward['src'], subject='sample', subjects_dir=subjects_dir, + clim=dict(kind='value', pos_lims=lims), mode='stat_map', + initial_time=0.1, verbose=True) ############################################################################### -# We can also visualize the activity on a "glass brain" (shown here with -# absolute values): +# Now let's plot this on a glass brain, which will automatically transform the +# data to MNI Talairach space: + +# sphinx_gallery_thumbnail_number = 4 -clim = dict(kind='value', lims=[0.3, 0.6, 0.9]) -abs(stc).plot(src=forward['src'], subject='sample', subjects_dir=subjects_dir, - mode='glass_brain', clim=clim) +stc.plot( + src=forward['src'], subject='sample', subjects_dir=subjects_dir, + mode='glass_brain', clim=dict(kind='value', lims=lims), + initial_time=0.1, verbose=True) + +############################################################################### +# Finally let's get another view, this time plotting again a ``'stat_map'`` +# style but using volumetric morphing to get data to fsaverage space, +# which we can get by passing a :class:`mne.SourceMorph` as the ``src`` +# argument to `mne.VolSourceEstimate.plot`. To save a bit of speed when +# applying the morph, we will crop the STC: + +morph = mne.compute_source_morph( + forward['src'], 'sample', 'fsaverage', subjects_dir=subjects_dir, + zooms=7, verbose=True) +stc.copy().crop(0.05, 0.18).plot( + src=morph, subject='fsaverage', subjects_dir=subjects_dir, + mode='stat_map', clim=dict(kind='value', pos_lims=lims), + initial_time=0.1, verbose=True) diff --git a/mne/conftest.py b/mne/conftest.py index 3d1fe80d9f6..b48b705fce2 100644 --- a/mne/conftest.py +++ b/mne/conftest.py @@ -73,6 +73,7 @@ def pytest_configure(config): ignore:numpy.ufunc size changed:RuntimeWarning ignore:.*mne-realtime.*:DeprecationWarning ignore:.*imp.*:DeprecationWarning + always:.*get_data.* is deprecated in favor of.*:DeprecationWarning """ # noqa: E501 for warning_line in warning_lines.split('\n'): warning_line = warning_line.strip() diff --git a/mne/fixes.py b/mne/fixes.py index 3d397caa141..251f0138b0f 100644 --- a/mne/fixes.py +++ b/mne/fixes.py @@ -420,6 +420,14 @@ def minimum_phase(h): ############################################################################### # Misc utilities +# Deal with nibabel 2.5 img.get_data() deprecation +def _get_img_fdata(img): + try: + return img.get_fdata() + except AttributeError: + return img.get_data().astype(float) + + def _read_volume_info(fobj): """An implementation of nibabel.freesurfer.io._read_volume_info, since old versions of nibabel (<=2.1.0) don't have it. diff --git a/mne/morph.py b/mne/morph.py index f1d62310e22..b684c78c04f 100644 --- a/mne/morph.py +++ b/mne/morph.py @@ -10,15 +10,16 @@ import numpy as np from scipy import sparse +from .fixes import _get_img_fdata from .parallel import parallel_func from .source_estimate import (VolSourceEstimate, SourceEstimate, VolVectorSourceEstimate, VectorSourceEstimate, _get_ico_tris) -from .source_space import SourceSpaces +from .source_space import SourceSpaces, _ensure_src from .surface import read_morph_map, mesh_edges, read_surface, _compute_nearest from .utils import (logger, verbose, check_version, get_subjects_dir, warn as warn_, deprecated, fill_doc, _check_option, - BunchConst) + BunchConst, wrapped_stdout) from .externals.h5io import read_hdf5, write_hdf5 @@ -119,6 +120,7 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage', kind = 'surface' subject_from = _check_subject_from(subject_from, src.subject) else: + src = _ensure_src(src) src_data, kind = _get_src_data(src) subject_from = _check_subject_from(subject_from, src) if not isinstance(subject_to, str): @@ -141,7 +143,7 @@ def compute_source_morph(src, subject_from=None, subject_to='fsaverage', shape = affine = pre_affine = sdr_morph = None morph_mat = vertices_to = None if kind == 'volume': - _check_dep(nibabel='2.1.0', dipy=False) + _check_dep(nibabel='2.1.0', dipy='0.10.1') logger.info('volume source space inferred...') import nibabel as nib @@ -324,7 +326,7 @@ def __init__(self, subject_from, subject_to, kind, zooms, @verbose def apply(self, stc_from, output='stc', mri_resolution=False, - mri_space=False, verbose=None): + mri_space=None, verbose=None): """Morph source space data. Parameters @@ -339,9 +341,9 @@ def apply(self, stc_from, output='stc', mri_resolution=False, If True the image is saved in MRI resolution. Default False. WARNING: if you have many time points the file produced can be huge. The default is mri_resolution=False. - mri_space : bool + mri_space : bool | None Whether the image to world registration should be in mri space. The - default is mri_space=mri_resolution. + default (None) is mri_space=mri_resolution. %(verbose_meth)s Returns @@ -468,8 +470,7 @@ def _check_dep(nibabel='2.1.0', dipy='0.10.1'): ver)) -def _morphed_stc_as_volume(morph, stc, mri_resolution=False, mri_space=True, - output='nifti1'): +def _morphed_stc_as_volume(morph, stc, mri_resolution, mri_space, output): """Return volume source space as Nifti1Image and/or save to disk.""" if isinstance(stc, VolVectorSourceEstimate): stc = stc.magnitude() @@ -478,15 +479,7 @@ def _morphed_stc_as_volume(morph, stc, mri_resolution=False, mri_space=True, 'volumes') _check_dep(nibabel='2.1.0', dipy=False) - _check_option('output', output, ['nifti', 'nifti1', 'nifti2']) - if output in ('nifti', 'nifti1'): - from nibabel import (Nifti1Image as NiftiImage, - Nifti1Header as NiftiHeader) - else: - assert output == 'nifti2' - from nibabel import (Nifti2Image as NiftiImage, - Nifti2Header as NiftiHeader) - + NiftiImage, NiftiHeader = _triage_output(output) new_zooms = None # if full MRI resolution, compute zooms from shape and MRI zooms @@ -512,7 +505,7 @@ def _morphed_stc_as_volume(morph, stc, mri_resolution=False, mri_space=True, img = np.zeros(morph.shape + (stc.shape[1],)).reshape(-1, stc.shape[1]) img[stc.vertices, :] = stc.data - img = img.reshape(morph.shape + (-1,)) + img = img.reshape(morph.shape + (-1,), order='F') # match order='F' above # make nifti from data with warnings.catch_warnings(): # nibabel<->numpy warning @@ -523,7 +516,7 @@ def _morphed_stc_as_volume(morph, stc, mri_resolution=False, mri_space=True, if new_zooms is not None: from dipy.align.reslice import reslice new_zooms = new_zooms[:3] - img, affine = reslice(img.get_data(), + img, affine = reslice(_get_img_fdata(img), img.affine, # MRI to world registration zooms, # old voxel size in mm new_zooms) # new voxel size in mm @@ -572,13 +565,8 @@ def _get_src_data(src): return src_data, src_kind -def _interpolate_data(stc, morph, mri_resolution=True, mri_space=True, - output='nifti1'): - """Interpolate source estimate data to MRI.""" - _check_dep(nibabel='2.1.0', dipy=False) - if output not in ('nifti', 'nifti1', 'nifti2'): - raise ValueError("invalid output specifier %s. Must be 'nifti1' or" - " 'nifti2'" % output) +def _triage_output(output): + _check_option('output', output, ['nifti', 'nifti1', 'nifti2']) if output in ('nifti', 'nifti1'): from nibabel import (Nifti1Image as NiftiImage, Nifti1Header as NiftiHeader) @@ -586,6 +574,13 @@ def _interpolate_data(stc, morph, mri_resolution=True, mri_space=True, assert output == 'nifti2' from nibabel import (Nifti2Image as NiftiImage, Nifti2Header as NiftiHeader) + return NiftiImage, NiftiHeader + + +def _interpolate_data(stc, morph, mri_resolution, mri_space, output): + """Interpolate source estimate data to MRI.""" + _check_dep(nibabel='2.1.0', dipy=False) + NiftiImage, NiftiHeader = _triage_output(output) assert morph.kind == 'volume' voxel_size_defined = False @@ -676,7 +671,8 @@ def _interpolate_data(stc, morph, mri_resolution=True, mri_space=True, if voxel_size_defined: # reslice mri img, img_affine = reslice( - img.get_data(), img.affine, _get_zooms_orig(morph), voxel_size) + _get_img_fdata(img), img.affine, _get_zooms_orig(morph), + voxel_size) with warnings.catch_warnings(): # nibabel<->numpy warning img = NiftiImage(img, img_affine, header=header) @@ -689,7 +685,6 @@ def _interpolate_data(stc, morph, mri_resolution=True, mri_space=True, def _compute_morph_sdr(mri_from, mri_to, niter_affine=(100, 100, 10), niter_sdr=(5, 5, 3), zooms=(5., 5., 5.)): """Get a matrix that morphs data from one subject to another.""" - _check_dep(nibabel='2.1.0', dipy='0.10.1') import nibabel as nib with np.testing.suppress_warnings(): from dipy.align import imaffine, imwarp, metrics, transforms @@ -709,25 +704,25 @@ def _compute_morph_sdr(mri_from, mri_to, niter_affine=(100, 100, 10), # reslice mri_from mri_from_res, mri_from_res_affine = reslice( - mri_from.get_data(), mri_from.affine, mri_from.header.get_zooms()[:3], - zooms) + _get_img_fdata(mri_from), mri_from.affine, + mri_from.header.get_zooms()[:3], zooms) with warnings.catch_warnings(): # nibabel<->numpy warning mri_from = nib.Nifti1Image(mri_from_res, mri_from_res_affine) # reslice mri_to mri_to_res, mri_to_res_affine = reslice( - mri_to.get_data(), mri_to.affine, mri_to.header.get_zooms()[:3], + _get_img_fdata(mri_to), mri_to.affine, mri_to.header.get_zooms()[:3], zooms) with warnings.catch_warnings(): # nibabel<->numpy warning mri_to = nib.Nifti1Image(mri_to_res, mri_to_res_affine) affine = mri_to.affine - mri_to = np.array(mri_to.dataobj, float) # to ndarray + mri_to = _get_img_fdata(mri_to) # to ndarray mri_to /= mri_to.max() mri_from_affine = mri_from.affine # get mri_from to world transform - mri_from = np.array(mri_from.dataobj, float) # to ndarray + mri_from = _get_img_fdata(mri_from) # to ndarray mri_from /= mri_from.max() # normalize # compute center of mass @@ -742,26 +737,33 @@ def _compute_morph_sdr(mri_from, mri_to, niter_affine=(100, 100, 10), factors=[4, 2, 1]) # translation - translation = affreg.optimize( - mri_to, mri_from, transforms.TranslationTransform3D(), None, affine, - mri_from_affine, starting_affine=c_of_mass.affine) + logger.info('Optimizing translation:') + with wrapped_stdout(indent=' '): + translation = affreg.optimize( + mri_to, mri_from, transforms.TranslationTransform3D(), None, + affine, mri_from_affine, starting_affine=c_of_mass.affine) # rigid body transform (translation + rotation) - rigid = affreg.optimize( - mri_to, mri_from, transforms.RigidTransform3D(), None, - affine, mri_from_affine, starting_affine=translation.affine) + logger.info('Optimizing rigid-body:') + with wrapped_stdout(indent=' '): + rigid = affreg.optimize( + mri_to, mri_from, transforms.RigidTransform3D(), None, + affine, mri_from_affine, starting_affine=translation.affine) # affine transform (translation + rotation + scaling) - pre_affine = affreg.optimize( - mri_to, mri_from, transforms.AffineTransform3D(), None, - affine, mri_from_affine, starting_affine=rigid.affine) + logger.info('Optimizing full affine:') + with wrapped_stdout(indent=' '): + pre_affine = affreg.optimize( + mri_to, mri_from, transforms.AffineTransform3D(), None, + affine, mri_from_affine, starting_affine=rigid.affine) # compute mapping sdr = imwarp.SymmetricDiffeomorphicRegistration( metrics.CCMetric(3), list(niter_sdr)) - sdr_morph = sdr.optimize(mri_to, pre_affine.transform(mri_from)) + logger.info('Optimizing SDR:') + with wrapped_stdout(indent=' '): + sdr_morph = sdr.optimize(mri_to, pre_affine.transform(mri_from)) shape = tuple(sdr_morph.domain_shape) # should be tuple of int - logger.info('done.') return shape, zooms, affine, pre_affine, sdr_morph @@ -1142,21 +1144,25 @@ def _apply_morph_data(morph, stc_from): def _morph_one(stc_one): # prepare data to be morphed + # here we use mri_resolution=True, mri_space=True because + # we will slice afterward + assert stc_one.data.shape[1] == 1 img_to = _interpolate_data(stc_one, morph, mri_resolution=True, - mri_space=True) + mri_space=True, output='nifti1') # reslice to match morph img_to, img_to_affine = reslice( - img_to.get_data(), morph.affine, _get_zooms_orig(morph), + _get_img_fdata(img_to), morph.affine, _get_zooms_orig(morph), morph.zooms) # morph data - for vol in range(img_to.shape[3]): - img_to[:, :, :, vol] = morph.sdr_morph.transform( - morph.pre_affine.transform(img_to[:, :, :, vol])) + img_to[:, :, :, 0] = morph.sdr_morph.transform( + morph.pre_affine.transform(img_to[:, :, :, 0])) - # reshape to nvoxel x nvol - img_to = img_to.reshape(-1, img_to.shape[3]) + # reshape to nvoxel x nvol: + # in the MNE definition of volume source spaces, + # x varies fastest, then y, then z, so we need order='F' here + img_to = img_to.reshape(-1, order='F') return img_to # First get the vertices (vertices_to) you will need the values for @@ -1165,7 +1171,7 @@ def _morph_one(stc_one): stc_from.vertices, tmin=0., tstep=1.) img_to = _morph_one(stc_ones) - vertices_to = np.where(img_to.sum(axis=1) != 0)[0] + vertices_to = np.where(img_to)[0] data = np.empty((len(vertices_to), n_times)) data_from = np.reshape(stc_from.data, (stc_from.data.shape[0], -1)) # Loop over time points to save memory @@ -1173,7 +1179,7 @@ def _morph_one(stc_one): this_stc = VolSourceEstimate( data_from[:, k:k + 1], stc_from.vertices, tmin=0., tstep=1.) this_img_to = _morph_one(this_stc) - data[:, k] = this_img_to[vertices_to, 0] + data[:, k] = this_img_to[vertices_to] data.shape = (len(vertices_to),) + stc_from.data.shape[1:] else: assert morph.kind == 'surface' diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 5a1cf612702..fcd3628b0fe 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -1331,7 +1331,7 @@ def to_original_src(self, src_orig, subject_orig=None, """ if self.subject is None: raise ValueError('stc.subject must be set') - src_orig = _ensure_src(src_orig, kind='surf') + src_orig = _ensure_src(src_orig, kind='surface') subject_orig = _ensure_src_subject(src_orig, subject_orig) data_idx, vertices = _get_morph_src_reordering( self.vertices, src_orig, subject_orig, self.subject, subjects_dir) @@ -1709,12 +1709,15 @@ def _vertices_list(self): @copy_function_doc_to_method_doc(plot_volume_source_estimates) def plot(self, src, subject=None, subjects_dir=None, mode='stat_map', bg_img=None, colorbar=True, colormap='auto', clim='auto', - transparent='auto', show=True, verbose=None): + transparent='auto', show=True, initial_time=None, + initial_pos=None, verbose=None): data = self.magnitude() if self._data_ndim == 3 else self return plot_volume_source_estimates( data, src=src, subject=subject, subjects_dir=subjects_dir, mode=mode, bg_img=bg_img, colorbar=colorbar, colormap=colormap, - clim=clim, transparent=transparent, show=show, verbose=verbose) + clim=clim, transparent=transparent, show=show, + initial_time=initial_time, initial_pos=initial_pos, + verbose=verbose) def save_as_volume(self, fname, src, dest='mri', mri_resolution=False, format='nifti1'): @@ -1777,7 +1780,7 @@ def as_volume(self, src, dest='mri', mri_resolution=False, Returns ------- - img : instance Nifti1Image + img : instance of Nifti1Image The image object. Notes @@ -2151,7 +2154,7 @@ def plot_surface(self, src, subject=None, surface='inflated', hemi='lh', A instance of `surfer.Brain` from PySurfer. """ # extract surface source spaces - surf = _ensure_src(src, kind='surf') + surf = _ensure_src(src, kind='surface') # extract surface source estimate data = self.data[:surf[0]['nuse'] + surf[1]['nuse']] @@ -2450,7 +2453,7 @@ def spatial_inter_hemi_connectivity(src, dist, verbose=None): using geodesic distances. """ from scipy.spatial.distance import cdist - src = _ensure_src(src, kind='surf') + src = _ensure_src(src, kind='surface') conn = cdist(src[0]['rr'][src[0]['vertno']], src[1]['rr'][src[1]['vertno']]) conn = sparse.csr_matrix(conn <= dist, dtype=int) diff --git a/mne/source_space.py b/mne/source_space.py index e0706d6ddb9..538942a6acc 100644 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -22,6 +22,7 @@ write_float_matrix, write_int_matrix, write_coord_trans, start_file, end_file, write_id) from .bem import read_bem_surfaces, ConductorModel +from .fixes import _get_img_fdata from .surface import (read_surface, _create_surf_spacing, _get_ico_surface, _tessellate_sphere_surf, _get_surf_neighbors, _normalize_vectors, _triangle_neighbors, mesh_dist, @@ -35,7 +36,7 @@ from .transforms import (invert_transform, apply_trans, _print_coord_trans, combine_transforms, _get_trans, _coord_frame_name, Transform, _str_to_frame, - _ensure_trans, _read_fs_xfm) + _ensure_trans, _read_ras_mni_t) def _get_lut(): @@ -354,7 +355,7 @@ def export_volume(self, fname, include_surfaces=True, # get the inuse array if mri_resolution: # read the mri file used to generate volumes - aseg_data = nib.load(vs['mri_file']).get_data() + aseg_data = _get_img_fdata(nib.load(vs['mri_file'])) # get the voxel space shape shape3d = (vs['mri_height'], vs['mri_depth'], vs['mri_width']) @@ -1272,18 +1273,15 @@ def head_to_mni(pos, subject, mri_head_t, subjects_dir=None, @verbose def _read_talxfm(subject, subjects_dir, mode=None, verbose=None): - """Read MNI transform from FreeSurfer talairach.xfm file. + """Compute MNI transform from FreeSurfer talairach.xfm file. Adapted from freesurfer m-files. Altered to deal with Norig and Torig correctly. """ if mode is not None and mode not in ['nibabel', 'freesurfer']: raise ValueError('mode must be "nibabel" or "freesurfer"') - fname = op.join(subjects_dir, subject, 'mri', 'transforms', - 'talairach.xfm') - # Setup the RAS to MNI transform - ras_mni_t = Transform('ras', 'mni_tal', _read_fs_xfm(fname)[0]) + ras_mni_t = _read_ras_mni_t(subject, subjects_dir) # We want to get from Freesurfer surface RAS ('mri') to MNI ('mni_tal'). # This file only gives us RAS (non-zero origin) ('ras') to MNI ('mni_tal'). @@ -1829,7 +1827,7 @@ def _get_volume_label_mask(mri, volume_label, rr): # Read the segmentation data using nibabel mgz = nib.load(mri) - mgz_data = mgz.get_data() + mgz_data = _get_img_fdata(mgz) # Get the numeric index for this volume label lut = _get_lut() @@ -2326,8 +2324,9 @@ def _adjust_patch_info(s, verbose=None): @verbose -def _ensure_src(src, kind=None, verbose=None): +def _ensure_src(src, kind=None, extra='', verbose=None): """Ensure we have a source space.""" + msg = 'src must be a string or instance of SourceSpaces%s' % (extra,) if _check_path_like(src): src = str(src) if not op.isfile(src): @@ -2335,14 +2334,10 @@ def _ensure_src(src, kind=None, verbose=None): logger.info('Reading %s...' % src) src = read_source_spaces(src, verbose=False) if not isinstance(src, SourceSpaces): - raise ValueError('src must be a string or instance of SourceSpaces') - if kind is not None: - if kind == 'surf': - surf = [s for s in src if s['type'] == 'surf'] - if len(surf) != 2 or len(src) != 2: - raise ValueError('Source space must contain exactly two ' - 'surfaces.') - src = surf + raise ValueError('%s, got %s (type %s)' % (msg, src, type(src))) + if kind is not None and src.kind != kind: + raise ValueError('Source space must be %s type, got ' + '%s' % (kind, src.kind)) return src @@ -2516,7 +2511,7 @@ def get_volume_labels_from_aseg(mgz_fname, return_colors=False): import nibabel as nib # Read the mgz file using nibabel - mgz_data = nib.load(mgz_fname).get_data() + mgz_data = _get_img_fdata(nib.load(mgz_fname)) # Get the unique label names lut = _get_lut() diff --git a/mne/tests/test_morph.py b/mne/tests/test_morph.py index b210ef10b5a..ebb7ed28af1 100644 --- a/mne/tests/test_morph.py +++ b/mne/tests/test_morph.py @@ -234,14 +234,14 @@ def test_volume_source_morph(): stc_vol = read_source_estimate(fname_vol, 'sample') # check for invalid input type - with pytest.raises(TypeError, match='src must be an instance of'): + with pytest.raises(ValueError, match='src must be a string or instance'): compute_source_morph(src=42) # check for raising an error if neither # inverse_operator_vol['src'][0]['subject_his_id'] nor subject_from is set, # but attempting to perform a volume morph src = inverse_operator_vol['src'] - src[0]['subject_his_id'] = None + assert src[0]['subject_his_id'] is None # already None on disk (old!) with pytest.raises(ValueError, match='subject_from could not be inferred'): compute_source_morph(src=src, subjects_dir=subjects_dir) @@ -259,7 +259,8 @@ def test_volume_source_morph(): zooms = 20 kwargs = dict(zooms=zooms, niter_sdr=(1,), niter_affine=(1,)) source_morph_vol = compute_source_morph( - subjects_dir=subjects_dir, src=inverse_operator_vol['src'], **kwargs) + subjects_dir=subjects_dir, src=fname_inv_vol, subject_from='sample', + **kwargs) shape = (13,) * 3 # for the given zooms assert source_morph_vol.subject_from == 'sample' diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 51122dfeda3..496657a1cc7 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -17,7 +17,7 @@ spatial_src_connectivity, spatial_tris_connectivity, SourceSpaces, VolVectorSourceEstimate) from mne.datasets import testing -from mne.fixes import fft +from mne.fixes import fft, _get_img_fdata from mne.source_estimate import grade_to_tris, _get_vol_mask from mne.minimum_norm import (read_inverse_operator, apply_inverse, apply_inverse_epochs) @@ -159,7 +159,7 @@ def test_stc_as_volume(): assert isinstance(img, nib.Nifti1Image) assert img.shape[:3] == inverse_operator_vol['src'][0]['shape'][:3] - with pytest.raises(ValueError, match='invalid output'): + with pytest.raises(ValueError, match='Invalid value.*output.*'): stc_vol.as_volume(inverse_operator_vol['src'], format='42') @@ -948,7 +948,7 @@ def test_vol_mask(): data = (1 + np.arange(n_vertices))[:, np.newaxis] stc_tmp = VolSourceEstimate(data, vertices, tmin=0., tstep=1.) img = stc_tmp.as_volume(src, mri_resolution=False) - img_data = img.get_data()[:, :, :, 0].T + img_data = _get_img_fdata(img)[:, :, :, 0].T mask_nib = (img_data != 0) assert_array_equal(img_data[mask_nib], data[:, 0]) assert_array_equal(np.where(mask_nib.ravel())[0], src[0]['vertno']) diff --git a/mne/transforms.py b/mne/transforms.py index ba86ba90370..bf7952034ae 100644 --- a/mne/transforms.py +++ b/mne/transforms.py @@ -19,7 +19,8 @@ from .io.open import fiff_open from .io.tag import read_tag from .io.write import start_file, end_file, write_coord_trans -from .utils import check_fname, logger, verbose, _ensure_int, _check_path_like +from .utils import (check_fname, logger, verbose, _ensure_int, _validate_type, + _check_path_like, get_subjects_dir) # transformation from anterior/left/superior coordinate system to @@ -173,14 +174,17 @@ def _coord_frame_name(cframe): return _verbose_frames.get(int(cframe), 'unknown') -def _print_coord_trans(t, prefix='Coordinate transformation: '): - logger.info(prefix + '%s -> %s' - % (_coord_frame_name(t['from']), _coord_frame_name(t['to']))) +def _print_coord_trans(t, prefix='Coordinate transformation: ', units='m', + level='info'): + # Units gives the units of the transformation. This always prints in mm. + log_func = getattr(logger, level) + log_func(prefix + '%s -> %s' + % (_coord_frame_name(t['from']), _coord_frame_name(t['to']))) for ti, tt in enumerate(t['trans']): - scale = 1000. if ti != 3 else 1. + scale = 1000. if (ti != 3 and units != 'mm') else 1. text = ' mm' if ti != 3 else '' - logger.info(' % 8.6f % 8.6f % 8.6f %7.2f%s' % - (tt[0], tt[1], tt[2], scale * tt[3], text)) + log_func(' % 8.6f % 8.6f % 8.6f %7.2f%s' % + (tt[0], tt[1], tt[2], scale * tt[3], text)) def _find_trans(subject, subjects_dir=None): @@ -1293,6 +1297,19 @@ def _average_quats(quats, weights=None): return avg_quat +def _read_ras_mni_t(subject, subjects_dir=None): + subjects_dir = get_subjects_dir(subjects_dir=subjects_dir, + raise_error=True) + _validate_type(subject, 'str', 'subject') + fname = op.join(subjects_dir, subject, 'mri', 'transforms', + 'talairach.xfm') + if not op.isfile(fname): + raise FileNotFoundError( + 'FreeSurfer Talairach transformation file not found: %s' + % (fname,)) + return Transform('ras', 'mni_tal', _read_fs_xfm(fname)[0]) + + def _read_fs_xfm(fname): """Read a Freesurfer transform from a .xfm file.""" assert fname.endswith('.xfm') @@ -1305,6 +1322,7 @@ def _read_fs_xfm(fname): for li, line in enumerate(fid): if li == 0: kind = line.strip() + logger.debug('Found: %r' % (kind,)) if line[:len(comp)] == comp: # we have the right line, so don't read any more break diff --git a/mne/utils/__init__.py b/mne/utils/__init__.py index 33fdc6e8846..9f8c7e83090 100644 --- a/mne/utils/__init__.py +++ b/mne/utils/__init__.py @@ -24,7 +24,7 @@ from .fetching import _fetch_file, _url_to_local_path from ._logging import (verbose, logger, set_log_level, set_log_file, use_log_level, catch_logging, warn, filter_out_warnings, - ETSContext) + ETSContext, wrapped_stdout) from .misc import (run_subprocess, _pl, _clean_names, pformat, _file_like, _explain_exception, _get_argvalues, sizeof_fmt, running_subprocess, _DefaultEventParser) diff --git a/mne/utils/_logging.py b/mne/utils/_logging.py index b302ae9a0bb..3d8c3aa1899 100644 --- a/mne/utils/_logging.py +++ b/mne/utils/_logging.py @@ -4,6 +4,7 @@ # # License: BSD (3-clause) +import contextlib import inspect from io import StringIO import re @@ -356,3 +357,17 @@ def __exit__(self, type, value, traceback): # noqa: D105 value.code += ("\nThis can probably be solved by setting " "ETS_TOOLKIT=qt4. On bash, type\n\n $ export " "ETS_TOOLKIT=qt4\n\nand run the command again.") + + +@contextlib.contextmanager +def wrapped_stdout(indent=''): + """Wrap stdout writes to logger.info, with an optional indent prefix.""" + orig_stdout = sys.stdout + my_out = StringIO() + sys.stdout = my_out + try: + yield + finally: + sys.stdout = orig_stdout + for line in my_out.getvalue().split('\n'): + logger.info(indent + line) diff --git a/mne/utils/check.py b/mne/utils/check.py index ec06fd2b54c..02f4ae6e355 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -149,20 +149,22 @@ def _check_fname(fname, overwrite=False, must_exist=False): return str(fname) -def _check_subject(class_subject, input_subject, raise_error=True): +def _check_subject(class_subject, input_subject, raise_error=True, + kind='class subject attribute'): """Get subject name from class.""" if input_subject is not None: _validate_type(input_subject, 'str', "subject input") + if class_subject is not None and input_subject != class_subject: + raise ValueError('%s (%r) did not match input subject (%r)' + % (kind, class_subject, input_subject)) return input_subject elif class_subject is not None: _validate_type(class_subject, 'str', - "Either subject input or class subject attribute") + "Either subject input or %s" % (kind,)) return class_subject - else: - if raise_error is True: - raise ValueError('Neither subject input nor class subject ' - 'attribute was a string') - return None + elif raise_error is True: + raise ValueError('Neither subject input nor %s was a string' % (kind,)) + return None def _check_preload(inst, msg): diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py index 6f7ff0e9b55..417e4e41763 100644 --- a/mne/viz/_3d.py +++ b/mne/viz/_3d.py @@ -22,20 +22,20 @@ from scipy import linalg, sparse from ..defaults import DEFAULTS -from ..fixes import einsum, _crop_colorbar +from ..fixes import einsum, _crop_colorbar, _get_img_fdata from ..io import _loc_to_coil_trans from ..io.pick import pick_types, _picks_to_idx from ..io.constants import FIFF from ..io.meas_info import read_fiducials, create_info from ..source_space import _ensure_src, _create_surf_spacing, _check_spacing -from ..surface import (get_meg_helmet_surf, read_surface, +from ..surface import (get_meg_helmet_surf, read_surface, _DistanceQuery, transform_surface_to, _project_onto_surface, - mesh_edges, _reorder_ccw, - _complete_sphere_surf) + mesh_edges, _reorder_ccw, _complete_sphere_surf) from ..transforms import (read_trans, _find_trans, apply_trans, rot_to_quat, combine_transforms, _get_trans, _ensure_trans, - invert_transform, Transform) + invert_transform, Transform, + _read_ras_mni_t, _print_coord_trans) from ..utils import (get_subjects_dir, logger, _check_subject, verbose, warn, has_nibabel, check_version, fill_doc, _pl, _ensure_int, _validate_type, _check_option) @@ -45,6 +45,7 @@ read_bem_surfaces) +verbose_dec = verbose FIDUCIAL_ORDER = (FIFF.FIFFV_POINT_LPA, FIFF.FIFFV_POINT_NASION, FIFF.FIFFV_POINT_RPA) @@ -431,7 +432,7 @@ def _plot_mri_contours(mri_fname, surf_fnames, orientation='coronal', # Load the T1 data nim = nib.load(mri_fname) - data = nim.get_data() + data = _get_img_fdata(nim) try: affine = nim.affine except AttributeError: # old nibabel @@ -1752,19 +1753,35 @@ def _glass_brain_crosshairs(params, x, y, z): ax.axhline(b, color='0.75') +def _cut_coords_to_ijk(cut_coords, img): + ijk = apply_trans(linalg.inv(img.affine), cut_coords) + ijk = np.clip(np.round(ijk).astype(int), 0, np.array(img.shape[:3]) - 1) + return ijk + + +def _ijk_to_cut_coords(ijk, img): + return apply_trans(img.affine, ijk) + + @verbose def plot_volume_source_estimates(stc, src, subject=None, subjects_dir=None, mode='stat_map', bg_img=None, colorbar=True, colormap='auto', clim='auto', - transparent=None, show=True, verbose=None): + transparent=None, show=True, + initial_time=None, initial_pos=None, + verbose=None): """Plot Nutmeg style volumetric source estimates using nilearn. Parameters ---------- stc : VectorSourceEstimate The vector source estimate to plot. - src : instance of SourceSpaces - The source space. + src : instance of SourceSpaces | instance of SourceMorph + The source space. Can also be a SourceMorph to morph the STC to + a new subject (see Examples). + + .. versionchanged:: 0.18 + Support for :class:`~nibabel.spatialimages.SpatialImage`. subject : str | None The subject name corresponding to FreeSurfer environment variable SUBJECT. If None stc.subject will be used. If that @@ -1776,7 +1793,7 @@ def plot_volume_source_estimates(stc, src, subject=None, subjects_dir=None, The plotting mode to use. Either 'stat_map' (default) or 'glass_brain'. For "glass_brain", activation absolute values are displayed after being transformed to a standard MNI brain. - bg_img : Niimg-like object | None + bg_img : instance of SpatialImage | None The background image used in the nilearn plotting function. If None, it is the T1.mgz file that is found in the subjects_dir. Not used in "glass brain" plotting. @@ -1787,6 +1804,20 @@ def plot_volume_source_estimates(stc, src, subject=None, subjects_dir=None, %(transparent)s show : bool Show figures if True. Defaults to True. + initial_time : float | None + The initial time to plot. Can be None (default) to use the time point + with the maximal absolute value activation across all voxels + or the ``initial_pos`` voxel (if ``initial_pos is None`` or not, + respectively). + + .. versionadded:: 0.19 + initial_pos : ndarray, shape (3,) | None + The initial position to use (in m). Can be None (default) to use the + voxel with the maximum absolute value activation across all time points + or at ``initial_time`` (if ``initial_time is None`` or not, + respectively). + + .. versionadded:: 0.19 %(verbose)s Notes @@ -1797,88 +1828,118 @@ def plot_volume_source_estimates(stc, src, subject=None, subjects_dir=None, The left and right arrow keys can be used to navigate in time. To move in time by larger steps, use shift+left and shift+right. + In ``'glass_brain'`` mode, values are transformed to the standard MNI + brain using the FreeSurfer Talairach transformation + ``$SUBJECTS_DIR/$SUBJECT/mri/transforms/talairach.xfm``. + .. versionadded:: 0.17 - """ + + .. versionchanged:: 0.19 + MRI volumes are automatically transformed to MNI space in + ``'glass_brain'`` mode. + + Examples + -------- + Passing a :class:`mne.SourceMorph` as the ``src`` + parameter can be useful for plotting in a different subject's space + (here, a ``'sample'`` STC in ``'fsaverage'``'s space):: + + >>> morph = mne.compute_source_morph(src_sample, subject_to='fsaverage') # doctest: +SKIP + >>> fig = stc_vol_sample.plot(morph) # doctest: +SKIP + + """ # noqa: E501 from matplotlib import pyplot as plt, colors from matplotlib.cbook import mplDeprecation import nibabel as nib from ..source_estimate import VolSourceEstimate + from ..morph import SourceMorph if not check_version('nilearn', '0.4'): raise RuntimeError('This function requires nilearn >= 0.4') from nilearn.plotting import plot_stat_map, plot_glass_brain - from nilearn.image import index_img, resample_to_img - - if mode == 'stat_map': - plot_func = plot_stat_map - elif mode == 'glass_brain': - plot_func = plot_glass_brain + from nilearn.image import index_img + + _check_option('mode', mode, ('stat_map', 'glass_brain')) + plot_func = dict(stat_map=plot_stat_map, + glass_brain=plot_glass_brain)[mode] + _validate_type(stc, VolSourceEstimate, 'stc') + if isinstance(src, SourceMorph): + img = src.apply(stc, 'nifti1', mri_resolution=False, mri_space=False) + stc = src.apply(stc, mri_resolution=False, mri_space=False) + kind, src_subject = 'morph.subject_to', src.subject_to else: - raise ValueError('Plotting function must be one of' - ' stat_map | glas_brain. Got %s' % mode) - - if not isinstance(stc, VolSourceEstimate): - raise ValueError('Only VolSourceEstimate objects are supported.' - 'Got %s' % type(stc)) - - def _cut_coords_to_idx(cut_coords, img_idx): + src = _ensure_src(src, kind='volume', extra=' or SourceMorph') + img = stc.as_volume(src, mri_resolution=False) + kind, src_subject = 'src subject', src[0].get('subject_his_id', None) + _print_coord_trans(Transform('mri_voxel', 'ras', img.affine), + prefix='Image affine ', units='mm', level='debug') + subject = _check_subject(src_subject, subject, True, kind=kind) + stc_ijk = np.array( + np.unravel_index(stc.vertices, img.shape[:3], order='F')).T + assert stc_ijk.shape == (len(stc.vertices), 3) + del src, kind + + # XXX this assumes zooms are uniform, should probably mult by zooms... + dist_to_verts = _DistanceQuery(stc_ijk, allow_kdtree=True) + + def _cut_coords_to_idx(cut_coords, img): """Convert voxel coordinates to index in stc.data.""" - # XXX: check lines below - cut_coords = apply_trans(linalg.inv(img.affine), cut_coords) - cut_coords = np.array([int(round(c)) for c in cut_coords]) - - # the affine transformation can sometimes lead to corner - # cases near the edges? - shape = img_idx.shape - cut_coords = np.clip(cut_coords, 0, np.array(shape) - 1) - loc_idx = np.ravel_multi_index( - cut_coords, shape, order='F') - dist_vertices = [abs(v - loc_idx) for v in stc.vertices] - nearest_idx = int(round(np.argmin(dist_vertices))) - if dist_vertices[nearest_idx] == 0: - return nearest_idx - else: - return None - - def _get_cut_coords_stat_map(event, params): + ijk = _cut_coords_to_ijk(cut_coords, img) + del cut_coords + logger.debug(' Affine remapped cut coords to [%d, %d, %d] idx' + % tuple(ijk)) + dist, loc_idx = dist_to_verts.query(ijk[np.newaxis]) + dist, loc_idx = dist[0], loc_idx[0] + logger.debug(' Using vertex %d at a distance of %d voxels' + % (stc.vertices[loc_idx], dist)) + return loc_idx + + ax_name = dict(x='X (saggital)', y='Y (coronal)', z='Z (axial)') + + def _click_to_cut_coords(event, params): """Get voxel coordinates from mouse click.""" if event.inaxes is params['ax_x']: - cut_coords = (params['ax_z'].lines[0].get_xdata()[0], - event.xdata, event.ydata) + ax = 'x' + x = params['ax_z'].lines[0].get_xdata()[0] + y, z = event.xdata, event.ydata elif event.inaxes is params['ax_y']: - cut_coords = (event.xdata, - params['ax_x'].lines[0].get_xdata()[0], - event.ydata) + ax = 'y' + y = params['ax_x'].lines[0].get_xdata()[0] + x, z = event.xdata, event.ydata + elif event.inaxes is params['ax_z']: + ax = 'z' + x, y = event.xdata, event.ydata + z = params['ax_x'].lines[1].get_ydata()[0] else: - cut_coords = (event.xdata, event.ydata, - params['ax_x'].lines[1].get_ydata()[0]) - return cut_coords + logger.debug(' Click outside axes') + return None + cut_coords = np.array((x, y, z)) + logger.debug('') + + if params['mode'] == 'glass_brain': # find idx for MIP + # Figure out what XYZ in world coordinates is in our voxel data + codes = ''.join(nib.aff2axcodes(params['img_idx'].affine)) + assert len(codes) == 3 + # We don't care about directionality, just which is which dim + codes = codes.replace('L', 'R').replace('P', 'A').replace('I', 'S') + idx = codes.index(dict(x='R', y='A', z='S')[ax]) + img_data = np.abs(_get_img_fdata(params['img_idx'])) + ijk = _cut_coords_to_ijk(cut_coords, params['img_idx']) + if idx == 0: + ijk[0] = np.argmax(img_data[:, ijk[1], ijk[2]]) + logger.debug(' MIP: i = %d idx' % (ijk[0],)) + elif idx == 1: + ijk[1] = np.argmax(img_data[ijk[0], :, ijk[2]]) + logger.debug(' MIP: j = %d idx' % (ijk[1],)) + else: + ijk[2] = np.argmax(img_data[ijk[0], ijk[1], :]) + logger.debug(' MIP: k = %d idx' % (ijk[2],)) + cut_coords = _ijk_to_cut_coords(ijk, params['img_idx']) - def _get_cut_coords_glass_brain(event, params): - """Get voxel coordinates with max intensity projection.""" - img_data = np.abs(params['img_idx_resampled'].get_data()) - shape = img_data.shape - if event.inaxes is params['ax_x']: - y, z = int(round(event.xdata)), int(round(event.ydata)) - x = np.argmax(img_data[:, y + shape[1] // 2, z + shape[2] // 2]) - x -= shape[0] // 2 - elif event.inaxes is params['ax_y']: - x, z = int(round(event.xdata)), int(round(event.ydata)) - y = np.argmax(img_data[x + shape[0] // 2, :, z + shape[2] // 2]) - y -= shape[1] // 2 - else: - x, y = int(round(event.xdata)), int(round(event.ydata)) - z = np.argmax(img_data[x + shape[0] // 2, y + shape[1] // 2, :]) - z -= shape[2] // 2 - return (x, y, z) - - def _resample(event, params): - """Precompute the resampling as the mouse leaves the time axis.""" - if event.inaxes is params['ax_time'] and mode == 'glass_brain': - img_resampled = resample_to_img(params['img_idx'], - params['bg_img']) - params.update({'img_idx_resampled': img_resampled}) + logger.debug(' Cut coords for %s: (%0.1f, %0.1f, %0.1f) mm' + % ((ax_name[ax],) + tuple(cut_coords))) + return cut_coords def _press(event, params): """Manage keypress on the plot.""" @@ -1892,19 +1953,21 @@ def _press(event, params): idx = min(params['stc'].shape[1] - 1, idx + 2) elif event.key == 'shift+right': idx = min(params['stc'].shape[1] - 1, idx + 10) - params['lx'].set_xdata(idx / params['stc'].sfreq + - params['stc'].tmin) _update_timeslice(idx, params) params['fig'].canvas.draw() def _update_timeslice(idx, params): - cut_coords = (0, 0, 0) + params['lx'].set_xdata(idx / params['stc'].sfreq + + params['stc'].tmin) ax_x, ax_y, ax_z = params['ax_x'], params['ax_y'], params['ax_z'] plot_map_callback = params['plot_func'] - if mode == 'stat_map': - cut_coords = (ax_y.lines[0].get_xdata()[0], - ax_x.lines[0].get_xdata()[0], - ax_x.lines[1].get_ydata()[0]) + # Crosshairs are the first thing plotted in stat_map, and the last + # in glass_brain + idxs = [0, 0, 1] if mode == 'stat_map' else [-2, -2, -1] + cut_coords = ( + ax_y.lines[idxs[0]].get_xdata()[0], + ax_x.lines[idxs[1]].get_xdata()[0], + ax_x.lines[idxs[2]].get_ydata()[0]) ax_x.clear() ax_y.clear() ax_z.clear() @@ -1912,69 +1975,123 @@ def _update_timeslice(idx, params): params.update({'title': 'Activation (t=%.3f s.)' % params['stc'].times[idx]}) plot_map_callback( - params['img_idx'], title='', - cut_coords=cut_coords) + params['img_idx'], title='', cut_coords=cut_coords) - def _onclick(event, params): + @verbose_dec + def _onclick(event, params, verbose=None): """Manage clicks on the plot.""" ax_x, ax_y, ax_z = params['ax_x'], params['ax_y'], params['ax_z'] plot_map_callback = params['plot_func'] if event.inaxes is params['ax_time']: - idx = params['stc'].time_as_index(event.xdata)[0] - params['lx'].set_xdata(event.xdata) + idx = params['stc'].time_as_index( + event.xdata, use_rounding=True)[0] _update_timeslice(idx, params) - if event.inaxes in [ax_x, ax_y, ax_z]: - if mode == 'stat_map': - cut_coords = _get_cut_coords_stat_map(event, params) - elif mode == 'glass_brain': - cut_coords = _get_cut_coords_glass_brain(event, params) - - ax_x.clear() - ax_y.clear() - ax_z.clear() - plot_map_callback(params['img_idx'], title='', - cut_coords=cut_coords) - loc_idx = _cut_coords_to_idx(cut_coords, params['img_idx']) - if loc_idx is not None: - ax_time.lines[0].set_ydata(stc.data[loc_idx].T) - else: - ax_time.lines[0].set_ydata(0.) - params['fig'].canvas.draw() + cut_coords = _click_to_cut_coords(event, params) + if cut_coords is None: + return # not in any axes - if bg_img is None: - subjects_dir = get_subjects_dir(subjects_dir=subjects_dir, - raise_error=True) - subject = _check_subject(stc.subject, subject, True) - t1_fname = op.join(subjects_dir, subject, 'mri', 'T1.mgz') - bg_img = nib.load(t1_fname) + ax_x.clear() + ax_y.clear() + ax_z.clear() + plot_map_callback(params['img_idx'], title='', + cut_coords=cut_coords) + loc_idx = _cut_coords_to_idx(cut_coords, params['img_idx']) + ydata = stc.data[loc_idx] + if loc_idx is not None: + ax_time.lines[0].set_ydata(ydata) + else: + ax_time.lines[0].set_ydata(0.) + params['fig'].canvas.draw() - bg_img_param = bg_img if mode == 'glass_brain': - bg_img_param = None - - img = stc.as_volume(src, mri_resolution=False) - - loc_idx, idx = np.unravel_index(np.abs(stc.data).argmax(), - stc.data.shape) - img_idx = index_img(img, idx) - x, y, z = np.unravel_index(stc.vertices[loc_idx], img_idx.shape, - order='F') - cut_coords = apply_trans(img.affine, (x, y, z)) + subject = _check_subject(stc.subject, subject, True) + ras_mni_t = _read_ras_mni_t(subject, subjects_dir) + if not np.allclose(ras_mni_t['trans'], np.eye(4)): + _print_coord_trans( + ras_mni_t, prefix='Transforming subject ', units='mm') + logger.info('') + # To get from voxel coords to world coords (i.e., define affine) + # we would apply img.affine, then also apply ras_mni_t, which + # transforms from the subject's RAS to MNI RAS. So we left-multiply + # these. + img = nib.Nifti1Image( + img.dataobj, np.dot(ras_mni_t['trans'], img.affine)) + bg_img = None # not used + else: # stat_map + if bg_img is None: + subject = _check_subject(stc.subject, subject, True) + subjects_dir = get_subjects_dir(subjects_dir=subjects_dir, + raise_error=True) + t1_fname = op.join(subjects_dir, subject, 'mri', 'T1.mgz') + bg_img = nib.load(t1_fname) + + if initial_time is None: + time_sl = slice(0, None) + else: + initial_time = float(initial_time) + logger.info('Fixing initial time: %s sec' % (initial_time,)) + initial_time = np.argmin(np.abs(stc.times - initial_time)) + time_sl = slice(initial_time, initial_time + 1) + if initial_pos is None: # find max pos and (maybe) time + loc_idx, time_idx = np.unravel_index( + np.abs(stc.data[:, time_sl]).argmax(), stc.data[:, time_sl].shape) + time_idx += time_sl.start + else: # position specified + initial_pos = np.array(initial_pos, float) + if initial_pos.shape != (3,): + raise ValueError('initial_pos must be float ndarray with shape ' + '(3,), got shape %s' % (initial_pos.shape,)) + initial_pos *= 1000 + logger.info('Fixing initial position: %s mm' + % (initial_pos.tolist(),)) + loc_idx = _cut_coords_to_idx(initial_pos, img) + if initial_time is not None: # time also specified + time_idx = time_sl.start + else: # find the max + time_idx = np.argmax(np.abs(stc.data[loc_idx])) + img_idx = index_img(img, time_idx) + assert img_idx.shape == img.shape[:3] + del initial_time, initial_pos + ijk = stc_ijk[loc_idx] + cut_coords = _ijk_to_cut_coords(ijk, img_idx) + np.testing.assert_allclose(_cut_coords_to_ijk(cut_coords, img_idx), ijk) + logger.info('Showing: t = %0.3f s, (%0.1f, %0.1f, %0.1f) mm, ' + '[%d, %d, %d] vox, %d vertex' + % ((stc.times[time_idx],) + tuple(cut_coords) + tuple(ijk) + + (stc.vertices[loc_idx],))) + del ijk # Plot initial figure fig, (axes, ax_time) = plt.subplots(2) - ax_time.plot(stc.times, stc.data[loc_idx].T, color='k') + axes.set(xticks=[], yticks=[]) + marker = 'o' if len(stc.times) == 1 else None + ydata = stc.data[loc_idx] + ax_time.plot(stc.times, ydata, color='k', marker=marker) if len(stc.times) > 1: ax_time.set(xlim=stc.times[[0, -1]]) ax_time.set(xlabel='Time (s)', ylabel='Activation') - lx = ax_time.axvline(stc.times[idx], color='g') - axes.set(xticks=[], yticks=[]) + lx = ax_time.axvline(stc.times[time_idx], color='g') fig.tight_layout() allow_pos_lims = (mode != 'glass_brain') - colormap, scale_pts, diverging, _, _ = _limits_to_control_points( + colormap, scale_pts, diverging, _, ticks = _limits_to_control_points( clim, stc.data, colormap, transparent, allow_pos_lims, linearize=True) + ylim = [min((scale_pts[0], ydata.min())), + max((scale_pts[-1], ydata.max()))] + ylim = np.array(ylim) + np.array([-1, 1]) * 0.05 * np.diff(ylim)[0] + dup_neg = False + if stc.data.min() < 0: + ax_time.axhline(0., color='0.5', ls='-', lw=0.5, zorder=2) + dup_neg = not diverging # glass brain with signed data + yticks = list(ticks) + if dup_neg: + yticks += [0] + list(-np.array(ticks)) + yticks = np.unique(yticks) + ax_time.set(yticks=yticks) + ax_time.set(ylim=ylim) + del yticks + if not diverging: # set eq above iff one-sided # there is a bug in nilearn where this messes w/transparency # Need to double the colormap @@ -2004,7 +2121,7 @@ def _onclick(event, params): plot_kwargs = dict( threshold=None, axes=axes, resampling_interpolation='nearest', vmax=vmax, figure=fig, - colorbar=colorbar, bg_img=bg_img_param, cmap=colormap, black_bg=True, + colorbar=colorbar, bg_img=bg_img, cmap=colormap, black_bg=True, symmetric_cbar=True) def plot_and_correct(*args, **kwargs): @@ -2015,6 +2132,7 @@ def plot_and_correct(*args, **kwargs): warnings.simplefilter('ignore', mplDeprecation) params['fig_anat'] = partial( plot_func, **plot_kwargs)(*args, **kwargs) + params['fig_anat']._cbar.outline.set_visible(False) for key in 'xyz': params.update({'ax_' + key: params['fig_anat'].axes[key].ax}) # Fix nilearn bug w/cbar background being white @@ -2023,26 +2141,22 @@ def plot_and_correct(*args, **kwargs): # adjust one-sided colorbars if not diverging: _crop_colorbar(params['fig_anat']._cbar, *scale_pts[[0, -1]]) + params['fig_anat']._cbar.set_ticks(params['cbar_ticks']) if mode == 'glass_brain': _glass_brain_crosshairs(params, *kwargs['cut_coords']) params = dict(stc=stc, ax_time=ax_time, plot_func=plot_and_correct, - img_idx=img_idx, fig=fig, bg_img=bg_img, lx=lx) + img_idx=img_idx, fig=fig, lx=lx, mode=mode, cbar_ticks=ticks) plot_and_correct(stat_map_img=params['img_idx'], title='', cut_coords=cut_coords) - if mode == 'glass_brain': - params.update(img_idx_resampled=resample_to_img( - params['img_idx'], params['bg_img'])) if show: plt.show() fig.canvas.mpl_connect('button_press_event', - partial(_onclick, params=params)) + partial(_onclick, params=params, verbose=verbose)) fig.canvas.mpl_connect('key_press_event', partial(_press, params=params)) - fig.canvas.mpl_connect('axes_leave_event', - partial(_resample, params=params)) return fig @@ -2670,7 +2784,7 @@ def _plot_dipole_mri_orthoview(dipole, trans, subject, subjects_dir=None, coord_frame=coord_frame) zooms = t1.header.get_zooms() - data = t1.get_data() + data = _get_img_fdata(t1) dims = len(data) # Symmetric size assumed. dd = dims / 2. dd *= t1.header.get_zooms()[0] diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 1cbaa41fd8e..b940810fbe6 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -19,7 +19,8 @@ read_trans, read_dipole, SourceEstimate, VectorSourceEstimate, VolSourceEstimate, make_sphere_model, use_coil_def, setup_volume_source_space, read_forward_solution, - VolVectorSourceEstimate, convert_forward_solution) + VolVectorSourceEstimate, convert_forward_solution, + compute_source_morph) from mne.io import read_raw_ctf, read_raw_bti, read_raw_kit, read_info from mne._digitization._utils import write_dig from mne.io.pick import pick_info @@ -475,41 +476,81 @@ def test_snapshot_brain_montage(renderer): @testing.requires_testing_data @requires_nibabel() @requires_version('nilearn', '0.4') -def test_plot_volume_source_estimates(): +@pytest.mark.parametrize('mode, stype, init_t, want_t, init_p, want_p', [ + ('glass_brain', 's', None, 2, None, (-30.9, 18.4, 56.7)), + ('stat_map', 'vec', 1, 1, None, (15.7, 16.0, -6.3)), + ('glass_brain', 'vec', None, 1, (10, -10, 20), (6.6, -9.0, 19.9)), + ('stat_map', 's', 1, 1, (-10, 5, 10), (-12.3, 2.0, 7.7))]) +def test_plot_volume_source_estimates(mode, stype, init_t, want_t, + init_p, want_p): """Test interactive plotting of volume source estimates.""" forward = read_forward_solution(fwd_fname) sample_src = forward['src'] + if init_p is not None: + init_p = np.array(init_p) / 1000. vertices = [s['vertno'] for s in sample_src] n_verts = sum(len(v) for v in vertices) n_time = 2 data = np.random.RandomState(0).rand(n_verts, n_time) - vol_stc = VolSourceEstimate(data, vertices, 1, 1) - - vol_vec_stc = VolVectorSourceEstimate( - np.tile(vol_stc.data[:, np.newaxis], (1, 3, 1)), vol_stc.vertices, - 0, 1) - for mode, stc in zip(['glass_brain', 'stat_map'], (vol_stc, vol_vec_stc)): - with pytest.warns(None): # sometimes get scalars/index warning - fig = stc.plot(sample_src, subject='sample', - subjects_dir=subjects_dir, - mode=mode) - # [ax_time, ax_y, ax_x, ax_z] - for ax_idx in [0, 2, 3, 4]: - _fake_click(fig, fig.axes[ax_idx], (0.3, 0.5)) - fig.canvas.key_press_event('left') - fig.canvas.key_press_event('shift+right') - - with pytest.raises(ValueError, match='must be one of'): - vol_stc.plot(sample_src, 'sample', subjects_dir, mode='abcd') + + if stype == 'vec': + stc = VolVectorSourceEstimate( + np.tile(data[:, np.newaxis], (1, 3, 1)), vertices, 1, 1) + else: + assert stype == 's' + stc = VolSourceEstimate(data, vertices, 1, 1) + with pytest.warns(None): # sometimes get scalars/index warning + with catch_logging() as log: + fig = stc.plot( + sample_src, subject='sample', subjects_dir=subjects_dir, + mode=mode, initial_time=init_t, initial_pos=init_p, + verbose=True) + log = log.getvalue() + want_str = 't = %0.3f s' % want_t + assert want_str in log, (want_str, init_t) + want_str = '(%0.1f, %0.1f, %0.1f) mm' % want_p + assert want_str in log, (want_str, init_p) + for ax_idx in [0, 2, 3, 4]: + _fake_click(fig, fig.axes[ax_idx], (0.3, 0.5)) + fig.canvas.key_press_event('left') + fig.canvas.key_press_event('shift+right') + + +@pytest.mark.slowtest # can be slow on OSX +@testing.requires_testing_data +@requires_nibabel() +@requires_version('nilearn', '0.4') +def test_plot_volume_source_estimates_morph(): + """Test interactive plotting of volume source estimates with morph.""" + forward = read_forward_solution(fwd_fname) + sample_src = forward['src'] + vertices = [s['vertno'] for s in sample_src] + n_verts = sum(len(v) for v in vertices) + n_time = 2 + data = np.random.RandomState(0).rand(n_verts, n_time) + stc = VolSourceEstimate(data, vertices, 1, 1) + morph = compute_source_morph(sample_src, 'sample', 'fsaverage', zooms=5, + subjects_dir=subjects_dir) + initial_pos = (-0.05, -0.01, -0.006) + with pytest.warns(None): # sometimes get scalars/index warning + with catch_logging() as log: + stc.plot(morph, subjects_dir=subjects_dir, mode='glass_brain', + initial_pos=initial_pos, verbose=True) + log = log.getvalue() + assert 't = 1.000 s' in log + assert '(-52.0, -8.0, -7.0) mm' in log + + with pytest.raises(ValueError, match='Allowed values are'): + stc.plot(sample_src, 'sample', subjects_dir, mode='abcd') vertices.append([]) surface_stc = SourceEstimate(data, vertices, 1, 1) - with pytest.raises(ValueError, match='Only Vol'): + with pytest.raises(TypeError, match='an instance of VolSourceEstimate'): plot_volume_source_estimates(surface_stc, sample_src, 'sample', subjects_dir) with pytest.raises(ValueError, match='Negative colormap limits'): - vol_stc.plot(sample_src, 'sample', subjects_dir, - clim=dict(lims=[-1, 2, 3], kind='value')) + stc.plot(sample_src, 'sample', subjects_dir, + clim=dict(lims=[-1, 2, 3], kind='value')) @testing.requires_testing_data