Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MRG, ENH: Add mri_resolution="sparse" #7888

Merged
merged 2 commits into from
Jun 11, 2020
Merged

MRG, ENH: Add mri_resolution="sparse" #7888

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4661319
ENH: Add mri_resolution="sparse"
larsoner Jun 10, 2020
f65d45f
FIX: Wording
larsoner Jun 11, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions doc/changes/latest.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`_
Expand Down Expand Up @@ -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`_
Expand Down
5 changes: 3 additions & 2 deletions examples/inverse/plot_mixed_source_space_inverse.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <[email protected]>
Expand Down Expand Up @@ -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)
Expand All @@ -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')

###############################################################################
Expand Down
2 changes: 1 addition & 1 deletion mne/forward/tests/test_make_forward.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
262 changes: 120 additions & 142 deletions mne/source_space.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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:
Expand Down Expand Up @@ -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
Expand Down
Loading