From 46613195b8d44c0eb679b1d6c5baa745cc65b5e0 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 10 Jun 2020 14:59:29 -0400 Subject: [PATCH 1/2] ENH: Add mri_resolution="sparse" --- doc/changes/latest.inc | 4 + .../plot_mixed_source_space_inverse.py | 5 +- mne/source_space.py | 262 ++++++++---------- mne/tests/test_source_space.py | 31 ++- 4 files changed, 152 insertions(+), 150 deletions(-) diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 7f5ded56550..721b0d612cd 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -73,6 +73,8 @@ Changelog - Add estimation method legend to :func:`mne.viz.plot_snr_estimate` by `Eric Larson`_ +- Add support to `mne.SourceSpaces.export_volume` for ``mri_resolution='sparse'`` to color only the nearest-neighbor voxels instead of entire regions by `Eric Larson`_ + - Add ``axes`` argument to :func:`mne.viz.plot_evoked_white`, :meth:`mne.Evoked.plot_white`, and :func:`mne.viz.plot_snr_estimate` by `Eric Larson`_ - Add ``sources`` and ``detectors`` options for fNIRS use of :meth:`mne.viz.plot_alignment` allowing plotting of optode locations in addition to channel midpoint ``channels`` and ``path`` between fNIRS optodes by `Kyle Mathewson`_ @@ -100,6 +102,8 @@ Bug - Fix bug with :func:`mne.setup_volume_source_space` when ``volume_label`` was supplied where voxels slightly (in a worst case, about 37% times ``pos`` in distance) outside the voxel-grid-based bounds of regions were errantly included, by `Eric Larson`_ +- Fix bug with `mne.SourceSpaces.export_volume` with ``use_lut=False`` where no values were written by `Eric Larson`_ + - Fix bug with :func:`mne.preprocessing.annotate_movement` where bad data segments, specified in ``raw.annotations``, would be handled incorrectly by `Luke Bloy`_ - Fix bug with :func:`mne.compute_source_morph` when more than one volume source space was present (e.g., when using labels) where only the first label would be interpolated when ``mri_resolution=True`` by `Eric Larson`_ diff --git a/examples/inverse/plot_mixed_source_space_inverse.py b/examples/inverse/plot_mixed_source_space_inverse.py index 04ed50e1471..81a5867697b 100644 --- a/examples/inverse/plot_mixed_source_space_inverse.py +++ b/examples/inverse/plot_mixed_source_space_inverse.py @@ -3,7 +3,7 @@ Compute MNE inverse solution on evoked data in a mixed source space =================================================================== -Create a mixed source space and compute MNE inverse solution on an +Create a mixed source space and compute an MNE inverse solution on an evoked dataset. """ # Author: Annalisa Pascarella @@ -79,6 +79,8 @@ print('the src space contains %d spaces and %d points' % (len(src), n)) ############################################################################### +# Viewing the source space +# ------------------------ # We could write the mixed source space with:: # # >>> write_source_spaces(fname_mixed_src, src, overwrite=True) @@ -87,7 +89,6 @@ nii_fname = op.join(bem_dir, '%s-mixed-src.nii' % subject) src.export_volume(nii_fname, mri_resolution=True, overwrite=True) - plotting.plot_img(nii_fname, cmap='nipy_spectral') ############################################################################### diff --git a/mne/source_space.py b/mne/source_space.py index 82224d9107f..8bc3e611d4e 100644 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -75,10 +75,8 @@ def _get_lut(fname=None): return np.genfromtxt(fname, dtype=dtype) -def _get_lut_id(lut, label, use_lut): +def _get_lut_id(lut, label): """Convert a label to a LUT ID number.""" - if not use_lut: - return 1 assert isinstance(label, str) mask = (lut['name'] == label) assert mask.sum() == 1 @@ -329,9 +327,15 @@ def export_volume(self, fname, include_surfaces=True, 4x4 transformation matrix (like the ``--trans`` MNE-C option. Must be provided if source spaces are in head coordinates and include_surfaces and mri_resolution are True. - mri_resolution : bool + mri_resolution : bool | str If True, the image is saved in MRI resolution - (e.g. 256 x 256 x 256). + (e.g. 256 x 256 x 256), and each source region (surface or + segmentation volume) filled in completely. If "sparse", only a + single voxel in the high-resolution MRI is filled in for each + source point. + + .. versionchanged:: 0.21.0 + Support for "sparse" was added. use_lut : bool If True, assigns a numeric value to each source space that corresponds to a color on the freesurfer lookup table. @@ -346,6 +350,12 @@ def export_volume(self, fname, include_surfaces=True, This method requires nibabel. """ _check_fname(fname, overwrite) + _validate_type(mri_resolution, (bool, str), 'mri_resolution') + if isinstance(mri_resolution, str): + _check_option('mri_resolution', mri_resolution, ["sparse"], + extra='when mri_resolution is a string') + else: + mri_resolution = bool(mri_resolution) fname = str(fname) # import nibabel or raise error try: @@ -388,170 +398,138 @@ def export_volume(self, fname, include_surfaces=True, lut = _get_lut() # Setup a dictionary of source types - src_types = dict(volume=[], surface=[], discrete=[]) + src_types = dict(volume=[], surface_discrete=[]) # Populate dictionary of source types for src in self: # volume sources if src['type'] == 'vol': src_types['volume'].append(src) - # surface sources - elif src['type'] == 'surf': - src_types['surface'].append(src) - # discrete sources - elif src['type'] == 'discrete': - src_types['discrete'].append(src) - # raise an error if dealing with source type other than volume - # surface or discrete + # surface and discrete sources + elif src['type'] in ('surf', 'discrete'): + src_types['surface_discrete'].append(src) else: raise ValueError('Unrecognized source type: %s.' % src['type']) + # Raise error if there are no volume source spaces + if len(src_types['volume']) == 0: + raise ValueError('Source spaces must contain at least one volume.') + # Get shape, inuse array and interpolation matrix from volume sources - inuse = 0 + src = src_types['volume'][0] + aseg_data = None + if mri_resolution: + # read the mri file used to generate volumes + if mri_resolution is True: + aseg_data = _get_img_fdata(nib.load(src['mri_file'])) + # get the voxel space shape + shape3d = (src['mri_width'], src['mri_depth'], + src['mri_height']) + else: + # get the volume source space shape + # read the shape in reverse order + # (otherwise results are scrambled) + shape3d = src['shape'] + + # calculate affine transform for image (MRI_VOXEL to RAS) + if mri_resolution: + # MRI_VOXEL to MRI transform + transform = src['vox_mri_t'] + else: + # MRI_VOXEL to MRI transform + # NOTE: 'src' indicates downsampled version of MRI_VOXEL + transform = src['src_mri_t'] + + # Figure out how to get from our input source space to output voxels + fro_dst_t = invert_transform(transform) + dest = transform['to'] + if coords == 'head': + head_mri_t = _get_trans(trans, 'head', 'mri')[0] + fro_dst_t = combine_transforms(head_mri_t, fro_dst_t, 'head', dest) + else: + fro_dst_t = fro_dst_t + + # Fill in the volumes + img = np.zeros(shape3d) for ii, vs in enumerate(src_types['volume']): # read the lookup table value for segmented volume if 'seg_name' not in vs: raise ValueError('Volume sources should be segments, ' 'not the entire volume.') # find the color value for this volume - id_ = _get_lut_id(lut, vs['seg_name'], use_lut) - - if ii == 0: - # get the inuse array - if mri_resolution: - # read the mri file used to generate volumes - 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']) - else: - # get the volume source space shape - # read the shape in reverse order - # (otherwise results are scrambled) - shape3d = vs['shape'][2::-1] - if mri_resolution: + use_id = 1. + if mri_resolution is True or use_lut: + id_ = _get_lut_id(lut, vs['seg_name']) + if use_lut: + use_id = id_ + + if mri_resolution == 'sparse': + idx = apply_trans(fro_dst_t, vs['rr'][vs['vertno']]) + idx = tuple(idx.round().astype(int).T) + elif mri_resolution is True: # fill the represented vol # get the values for this volume - use = id_ * (aseg_data == id_).astype(int).ravel('F') + idx = (aseg_data == id_) else: - use = id_ * vs['inuse'] - inuse += use - - # Raise error if there are no volume source spaces - if np.array(inuse).ndim == 0: - raise ValueError('Source spaces must contain at least one volume.') - - # create 3d grid in the MRI_VOXEL coordinate frame - # len of inuse array should match shape regardless of mri_resolution - assert len(inuse) == np.prod(shape3d) - - # setup the image in 3d space - img = inuse.reshape(shape3d).T - - # include surface and/or discrete source spaces - if include_surfaces or include_discrete: - - # setup affine transform for source spaces - if mri_resolution: - # get the MRI to MRI_VOXEL transform - affine = invert_transform(vs['vox_mri_t']) + assert mri_resolution is False + idx = vs['inuse'].reshape(shape3d, order='F').astype(bool) + img[idx] = use_id + + # loop through the surface and discrete source spaces + + # get the surface names (assumes left, right order. may want + # to add these names during source space generation + for src in src_types['surface_discrete']: + val = 1 + if src['type'] == 'surf': + if not include_surfaces: + continue + if use_lut: + surf_name = { + FIFF.FIFFV_MNE_SURF_LEFT_HEMI: 'Left', + FIFF.FIFFV_MNE_SURF_RIGHT_HEMI: 'Right', + }[src['id']] + '-Cerebral-Cortex' + val = _get_lut_id(lut, surf_name) else: - # get the MRI to SOURCE (MRI_VOXEL) transform - affine = invert_transform(vs['src_mri_t']) - - # modify affine if in head coordinates - if coords == 'head': - - # read mri -> head transformation - mri_head_t = _get_trans(trans)[0] - - # get the HEAD to MRI transform - head_mri_t = invert_transform(mri_head_t) - - # combine transforms, from HEAD to MRI_VOXEL - affine = combine_transforms(head_mri_t, affine, - 'head', 'mri_voxel') - - # loop through the surface source spaces - if include_surfaces: - - # get the surface names (assumes left, right order. may want - # to add these names during source space generation - surf_names = ['Left-Cerebral-Cortex', 'Right-Cerebral-Cortex'] - - for i, surf in enumerate(src_types['surface']): - # convert vertex positions from their native space - # (either HEAD or MRI) to MRI_VOXEL space - srf_rr = apply_trans(affine['trans'], surf['rr']) - # convert to numeric indices - ix_orig, iy_orig, iz_orig = srf_rr.T.round().astype(int) - # clip indices outside of volume space - ix_clip = np.maximum(np.minimum(ix_orig, shape3d[2] - 1), - 0) - iy_clip = np.maximum(np.minimum(iy_orig, shape3d[1] - 1), - 0) - iz_clip = np.maximum(np.minimum(iz_orig, shape3d[0] - 1), - 0) - # compare original and clipped indices - n_diff = np.array((ix_orig != ix_clip, iy_orig != iy_clip, - iz_orig != iz_clip)).any(0).sum() - # generate use warnings for clipping - if n_diff > 0: - warn('%s surface vertices lay outside of volume space.' - ' Consider using a larger volume space.' % n_diff) - # get surface id or use default value - i = _get_lut_id(lut, surf_names[i], use_lut) - # update image to include surface voxels - img[ix_clip, iy_clip, iz_clip] = i - - # loop through discrete source spaces - if include_discrete: - for i, disc in enumerate(src_types['discrete']): - # convert vertex positions from their native space - # (either HEAD or MRI) to MRI_VOXEL space - disc_rr = apply_trans(affine['trans'], disc['rr']) - # convert to numeric indices - ix_orig, iy_orig, iz_orig = disc_rr.T.astype(int) - # clip indices outside of volume space - ix_clip = np.maximum(np.minimum(ix_orig, shape3d[2] - 1), - 0) - iy_clip = np.maximum(np.minimum(iy_orig, shape3d[1] - 1), - 0) - iz_clip = np.maximum(np.minimum(iz_orig, shape3d[0] - 1), - 0) - # compare original and clipped indices - n_diff = np.array((ix_orig != ix_clip, iy_orig != iy_clip, - iz_orig != iz_clip)).any(0).sum() - # generate use warnings for clipping - if n_diff > 0: - warn('%s discrete vertices lay outside of volume ' - 'space. Consider using a larger volume space.' - % n_diff) - # set default value - img[ix_clip, iy_clip, iz_clip] = 1 - if use_lut: - logger.info('Discrete sources do not have values on ' - 'the lookup table. Defaulting to 1.') + assert src['type'] == 'discrete' + if not include_discrete: + continue + if use_lut: + logger.info('Discrete sources do not have values on ' + 'the lookup table. Defaulting to 1.') + # convert vertex positions from their native space + # (either HEAD or MRI) to MRI_VOXEL space + if mri_resolution is True: + use_rr = src['rr'] + else: + assert mri_resolution is False or mri_resolution == 'sparse' + use_rr = src['rr'][src['vertno']] + srf_vox = apply_trans(fro_dst_t['trans'], use_rr) + # convert to numeric indices + ix_, iy_, iz_ = srf_vox.T.round().astype(int) + # clip indices outside of volume space + ix = np.clip(ix_, 0, shape3d[0] - 1), + iy = np.clip(iy_, 0, shape3d[1] - 1) + iz = np.clip(iz_, 0, shape3d[2] - 1) + # compare original and clipped indices + n_diff = ((ix_ != ix) | (iy_ != iy) | (iz_ != iz)).sum() + # generate use warnings for clipping + if n_diff > 0: + warn(f'{n_diff} {src["kind"]} vertices lay outside of volume ' + f'space. Consider using a larger volume space.') + # get surface id or use default value + # update image to include surface voxels + img[ix, iy, iz] = val - # calculate affine transform for image (MRI_VOXEL to RAS) - if mri_resolution: - # MRI_VOXEL to MRI transform - transform = vs['vox_mri_t'] - else: - # MRI_VOXEL to MRI transform - # NOTE: 'src' indicates downsampled version of MRI_VOXEL - transform = vs['src_mri_t'] if dest == 'mri': # combine with MRI to RAS transform - transform = combine_transforms(transform, vs['mri_ras_t'], - transform['from'], - vs['mri_ras_t']['to']) + transform = combine_transforms( + transform, vs['mri_ras_t'], + transform['from'], vs['mri_ras_t']['to']) # now setup the affine for volume image affine = transform['trans'].copy() # make sure affine converts from m to mm affine[:3] *= 1e3 - # save volume data - # setup image for file if fname.endswith(('.nii', '.nii.gz')): # save as nifit # setup the nifti header diff --git a/mne/tests/test_source_space.py b/mne/tests/test_source_space.py index 720514cb5ce..eafb56e5bd8 100644 --- a/mne/tests/test_source_space.py +++ b/mne/tests/test_source_space.py @@ -19,6 +19,7 @@ morph_source_spaces, SourceEstimate, make_sphere_model, head_to_mni, read_trans, compute_source_morph, read_bem_solution, read_freesurfer_lut) +from mne.fixes import _get_img_fdata from mne.utils import (requires_nibabel, run_subprocess, modified_env, requires_mne, run_tests_if_main, check_version) @@ -682,13 +683,15 @@ def test_read_volume_from_src(): @requires_nibabel() def test_combine_source_spaces(tmpdir): """Test combining source spaces.""" - rng = np.random.RandomState(0) + import nibabel as nib + rng = np.random.RandomState(2) aseg_fname = op.join(subjects_dir, 'sample', 'mri', 'aseg.mgz') - label_names = get_volume_labels_from_aseg(aseg_fname) - volume_labels = rng.choice(label_names, 2) + volume_labels = ['Brain-Stem', 'Right-Hippocampus'] # two fairly large - # get a surface source space (no need to test creation here) - srf = read_source_spaces(fname, patch_stats=False) + # create a sparse surface source space to ensure all get mapped + # when mri_resolution=False + srf = setup_source_space('sample', 'oct3', add_dist=False, + subjects_dir=subjects_dir) # setup 2 volume source spaces vol = setup_volume_source_space('sample', subjects_dir=subjects_dir, @@ -696,7 +699,7 @@ def test_combine_source_spaces(tmpdir): mri=aseg_fname, add_interpolator=False) # setup a discrete source space - rr = rng.randint(0, 20, (100, 3)) * 1e-3 + rr = rng.randint(0, 11, (20, 3)) * 5e-3 nn = np.zeros(rr.shape) nn[:, -1] = 1 pos = {'rr': rr, 'nn': nn} @@ -754,6 +757,22 @@ def test_combine_source_spaces(tmpdir): pytest.raises(ValueError, src_mixed_coord.export_volume, image_fname, verbose='error') + # now actually write it + fname_img = tmpdir.join('img.nii') + for mri_resolution in (False, 'sparse', True): + for src, up in ((vol, 705), + (srf + vol, 27272), + (disc + vol, 705)): + src.export_volume( + fname_img, use_lut=False, + mri_resolution=mri_resolution, overwrite=True) + img_data = _get_img_fdata(nib.load(str(fname_img))) + n_src = img_data.astype(bool).sum() + n_want = sum(s['nuse'] for s in src) + if mri_resolution is True: + n_want += up + assert n_src == n_want, src + @testing.requires_testing_data def test_morph_source_spaces(): From f65d45fefe197cefb6c8b379e2f477048405f762 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 11 Jun 2020 11:22:34 -0400 Subject: [PATCH 2/2] FIX: Wording --- mne/forward/tests/test_make_forward.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index f145e22b708..acab22df01d 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -336,7 +336,7 @@ def test_forward_mixed_source_space(tmpdir): # head coordinates and mri_resolution, but wrong trans file vox_mri_t = vol1[0]['vox_mri_t'] - with pytest.raises(ValueError, match='mri<->head, got mri_voxel->mri'): + with pytest.raises(ValueError, match='head<->mri, got mri_voxel->mri'): src_from_fwd.export_volume(fname_img, mri_resolution=True, trans=vox_mri_t)