Skip to content

Commit

Permalink
major coordinate frame reworking, almost there pending misc data upda…
Browse files Browse the repository at this point in the history
…te, unknown PR
  • Loading branch information
alexrockhill committed Jul 30, 2021
1 parent f0a6981 commit c1d78df
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 104 deletions.
28 changes: 13 additions & 15 deletions mne/_freesurfer.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@
from .io.constants import FIFF
from .io.meas_info import read_fiducials
from .transforms import (apply_trans, invert_transform, combine_transforms,
_ensure_trans, read_ras_mni_t, Transform, _get_trans)
_ensure_trans, read_ras_mni_t, Transform)
from .surface import read_surface, _read_mri_surface
from .utils import (verbose, _validate_type, _check_fname, _check_option,
get_subjects_dir, _require_version, logger)
get_subjects_dir, _require_version, logger, warn)


def _check_subject_dir(subject, subjects_dir):
Expand Down Expand Up @@ -445,8 +445,9 @@ def get_mri_head_trans(subject, subjects_dir, verbose=None):
return trans


def _coord_to_mri(info, idx, subject, subjects_dir, trans=None,
sensor=False, detector=False):
@verbose
def _coord_to_mri(info, idx, subject, subjects_dir,
sensor=False, detector=False, verbose=None):
"""Transform a channel to the "mri"/surface RAS coordinate frame."""
assert sensor + detector < 2 # can't be both
type_slice = slice(0, 3)
Expand All @@ -458,24 +459,21 @@ def _coord_to_mri(info, idx, subject, subjects_dir, trans=None,
# coordinate frame/location info
coord_frame = info['chs'][idx]['coord_frame']
if coord_frame == FIFF.FIFFV_COORD_UNKNOWN:
if trans is None: # need to apply user-supplied trans
raise RuntimeError('`trans` must not be `None` when '
'the coordinate frame is "unknown"')
head_mri_t = _get_trans(trans, 'head', 'mri')[0]
ch_coord = apply_trans(trans['trans'], ch_coord)
coord_frame = trans['to']
if verbose is None or verbose or verbose == 'warn':
warn('Got coordinate frame "unknown", assuming '
'"head" coordinates')
coord_frame = FIFF.FIFFV_COORD_HEAD
if coord_frame == FIFF.FIFFV_COORD_DEVICE:
ch_coord = apply_trans(info['dev_head_t']['trans'], ch_coord)
coord_frame = FIFF.FIFFV_COORD_HEAD # now it's in head
if coord_frame == FIFF.FIFFV_COORD_HEAD:
trans = get_mri_head_trans(subject, subjects_dir)
ch_coord = apply_trans(trans['trans'], ch_coord) * 1000
trans = invert_transform(get_mri_head_trans(subject, subjects_dir))
ch_coord = apply_trans(trans['trans'], ch_coord) * 1000 # mm -> m
coord_frame = FIFF.FIFFV_COORD_MRI # now in surface RAS
if coord_frame != FIFF.FIFFV_COORD_MRI:
raise RuntimeError(
'`inst` coordinate frame must be "unknown", "meg", '
f'"head" or "mri", got {coord_frame} for '
f'{info.ch_names[idx]}.')
'`inst` coordinate frame must be "meg", "head" or "mri", '
f'got {coord_frame} for {info.ch_names[idx]}')
return ch_coord


Expand Down
48 changes: 17 additions & 31 deletions mne/surface.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
from .fixes import (_serialize_volume_info, _get_read_geometry, jit,
prange, bincount, _get_img_fdata)
from .io.constants import FIFF
from .io.meas_info import Info
from .io.pick import pick_types
from .parallel import parallel_func
from .transforms import (transform_surface_to, _pol_to_cart, _cart_to_sph,
Expand Down Expand Up @@ -1752,19 +1753,19 @@ def _get_surface_RAS_volumes(base_image, subject_from, subjects_dir):
return base_image, fs_t1


def _warn_missing_chs(montage, dig_image, after_warp):
def _warn_missing_chs(info, dig_image, after_warp):
"""Warn that channels are missing."""
# ensure that each electrode contact was marked in at least one voxel
missing = set(np.arange(1, len(montage.ch_names) + 1)).difference(
missing = set(np.arange(1, len(info.ch_names) + 1)).difference(
set(np.unique(_get_img_fdata(dig_image))))
missing_ch = [montage.ch_names[idx - 1] for idx in missing]
missing_ch = [info.ch_names[idx - 1] for idx in missing]
if missing_ch:
warn('Channels ' + ', '.join(missing_ch) + ' were not assigned '
'voxels' + (' after applying SDR warp' if after_warp else ''))


@fill_doc
def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
def warp_montage_volume(info, base_image, reg_affine, sdr_morph,
subject_from, subject_to='fsaverage',
subjects_dir=None, thresh=0.95,
max_peak_dist=1, voxels_max=100, use_min=False):
Expand All @@ -1780,8 +1781,7 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
Parameters
----------
montage : mne.channels.montage.DigMontage
The montage object containing the channels.
%(info)s
base_image : str | pathlib.Path | nibabel.spatialimages.SpatialImage
Path to a volumetric scan (e.g. CT) of the subject. Can be in any
format readable by nibabel. Can also be a nibabel image object.
Expand Down Expand Up @@ -1810,8 +1810,6 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
Returns
-------
montage_warped : mne.channels.montage.DigMontage
The modified montage object containing the channels.
image_from : nibabel.spatialimages.SpatialImage
An image in Freesurfer surface RAS space with voxel values
corresponding to the index of the channel. The background
Expand All @@ -1822,11 +1820,10 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
"""
_require_version('nibabel', 'SDR morph', '2.1.0')
_require_version('dipy', 'SDR morph', '0.10.1')
from .channels import DigMontage
from ._freesurfer import _check_subject_dir
from ._freesurfer import _check_subject_dir, _coord_to_mri
import nibabel as nib

_validate_type(montage, DigMontage, 'montage')
_validate_type(info, Info, 'info')
_validate_type(base_image, nib.spatialimages.SpatialImage, 'base_image')
_validate_type(thresh, float, 'thresh')
if thresh < 0 or thresh >= 1:
Expand All @@ -1839,17 +1836,13 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
_check_subject_dir(subject_from, subjects_dir)
_check_subject_dir(subject_to, subjects_dir)

ch_coords = np.array(
[_coord_to_mri(info, idx) for idx in range(info['nchan'])])

# load image and make sure it's in surface RAS
image, fs_t1 = _get_surface_RAS_volumes(
base_image, subject_from, subjects_dir)

# get montage channel coordinates
ch_dict = montage.get_positions()
if ch_dict['coord_frame'] != 'mri':
raise RuntimeError('Coordinate frame not supported, expected '
'"mri", got ' + str(ch_dict['coord_frame']))
ch_coords = np.array(list(ch_dict['ch_pos'].values()))

# convert to freesurfer voxel space
ch_coords = apply_trans(
np.linalg.inv(fs_t1.header.get_vox2ras_tkr()), ch_coords * 1000)
Expand Down Expand Up @@ -1881,7 +1874,7 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,

# apply the mapping
image_from = nib.Nifti1Image(image_from, fs_t1.affine)
_warn_missing_chs(montage, image_from, after_warp=False)
_warn_missing_chs(info, image_from, after_warp=False)

template_brain = nib.load(
op.join(subjects_dir, subject_to, 'mri', 'brain.mgz'))
Expand All @@ -1890,7 +1883,7 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
image_from, template_brain, reg_affine, sdr_morph,
interpolation='nearest')

_warn_missing_chs(montage, image_to, after_warp=True)
_warn_missing_chs(info, image_to, after_warp=True)

# recover the contact positions as the center of mass
warped_data = _get_img_fdata(image_to)
Expand All @@ -1902,17 +1895,10 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
ch_coords = apply_trans(
fs_t1.header.get_vox2ras_tkr(), ch_coords) / 1000

# copy before modifying
montage_warped = montage.copy()

# modify montage to be returned
idx = 0
for dig_point in montage_warped.dig:
if dig_point['kind'] == FIFF.FIFFV_POINT_EEG:
assert dig_point['ident'] == idx + 1 # 1 indexed
dig_point['r'] = ch_coords[idx]
idx += 1
return montage_warped, image_from, image_to
# modify info in place
for ch, ch_coord in zip(info['chs'], ch_coords):
ch['loc'][:3] = ch_coords
return image_from, image_to


_VOXELS_MAX = 100 # define constant to avoid runtime issues
Expand Down
19 changes: 9 additions & 10 deletions mne/viz/_brain/_brain.py
Original file line number Diff line number Diff line change
Expand Up @@ -2502,18 +2502,17 @@ def add_foci(self, coords, coords_as_verts=False, map_surface=None,
opacity=alpha, resolution=resolution)
self._renderer.set_camera(**views_dicts[hemi][v])

@fill_doc
def add_sensors(self, inst, trans=None, picks=None, colors=None,
@verbose
def add_sensors(self, inst, picks=None, colors=None,
scales=None, alphas=None, meg_helmet=True,
project_eeg=False, fnirs_pairs='gray',
fnirs_sources='red', fnirs_detectors='brown'):
fnirs_sources='red', fnirs_detectors='brown', verbose=None):
"""Add mesh objects to represent sensor positions.
Parameters
----------
inst : mne.io.Raw | mne.Epochs | mne.Evoked
The MNE object with sensor locations and data.
%(trans)s
%(picks_all)s
colors : str | dict
The color of the sensors. If dict, channel names or types can be
Expand Down Expand Up @@ -2542,6 +2541,7 @@ def add_sensors(self, inst, trans=None, picks=None, colors=None,
fnirs_detectors : str | bool
The color to plot the fnirs detectors. If False, fnirs detectors
will not be plotted, if True, the color defaults to "brown".
%(verbose)s
Notes
-----
Expand All @@ -2552,14 +2552,13 @@ def add_sensors(self, inst, trans=None, picks=None, colors=None,
allow_empty=False)
info = inst.info.copy().pick_channels(
[inst.ch_names[idx] for idx in picks])
montage = inst.get_montage()
if project_eeg and pick_types(info.copy(), eeg=True).size > 0:
head_surf = _get_head_surface(
'seghead', self._subject_id, self._subjects_dir)
for idx in range(info['nchan']):
try:
ch_coord = _coord_to_mri(info, idx, self._subject_id,
self._subjects_dir, trans)
self._subjects_dir, verbose=verbose)
except RuntimeError as e:
raise RuntimeError(
f'{e} You may want to use `picks=("meg", "eeg")` '
Expand All @@ -2582,11 +2581,11 @@ def add_sensors(self, inst, trans=None, picks=None, colors=None,
center=tuple(ch_coord), color=color, scale=scale)
if ch_type in _FNIRS_CH_TYPES_SPLIT:
sensor_ch_coord = _coord_to_mri(
info, idx, self._subject_id, self._subjects_dir, trans,
sensor=True)
info, idx, self._subject_id, self._subjects_dir,
sensor=True, verbose=verbose)
detector_ch_coord = _coord_to_mri(
info, idx, self._subject_id, self._subjects_dir, trans,
detector=True)
info, idx, self._subject_id, self._subjects_dir,
detector=True, verbose=verbose)
if fnirs_pairs:
color = colorConverter.to_rgba(
'gray' if isinstance(fnirs_pairs, bool) else
Expand Down
47 changes: 11 additions & 36 deletions tutorials/clinical/10_ieeg_localize.py
Original file line number Diff line number Diff line change
Expand Up @@ -205,13 +205,6 @@ def plot_overlay(image, compare, title, thresh=None):
# load electrophysiology data with channel locations
raw = mne.io.read_raw(op.join(misc_path, 'seeg', 'sample_seeg_ieeg.fif'))

# MNE stores digitization montages in a coordinate frame called "head"
# thus even though the channel positions were found in the reference frame of
# Freesurfer surfaces when they were assigned to the ``raw`` object, they were
# transformed to "head" so we need to have a transform to put them back
subj_trans = mne.read_trans(op.join(misc_path, 'seeg',
'sample_seeg_trans.fif'))

# create symbolic link to share ``subjects_dir``
seeg_subject_dir = op.join(subjects_dir, 'sample_seeg')
if not op.exists(seeg_subject_dir):
Expand All @@ -220,19 +213,13 @@ def plot_overlay(image, compare, title, thresh=None):
###############################################################################
# Let's plot the electrode contact locations on the subject's brain.

# load the subject's brain
subject_brain = nib.load(
op.join(misc_path, 'seeg', 'sample_seeg', 'mri', 'brain.mgz'))

# plot the alignment
fig_kwargs = dict(size=(800, 600), bgcolor='w', scene=False)
renderer = mne.viz.backends.renderer.create_3d_figure(**fig_kwargs)
brain = mne.viz.Brain('sample_seeg', hemi='both', surf='pial',
cortex='low_contrast', alpha=0.2,
subjects_dir=subjects_dir, figure=renderer.figure)
brain_kwargs = dict(cortex='low_contrast', alpha=0.1,
subjects_dir=subjects_dir)
brain = mne.viz.Brain('sample_seeg', **brain_kwargs)
brain.add_sensors(raw)
view_kwargs = dict(azimuth=60, elevation=100, distance=300)
mne.viz.set_3d_view(renderer.figure, **view_kwargs)
brain.show_view(**view_kwargs)

# %%
# Warping to a Common Atlas
Expand All @@ -248,7 +235,9 @@ def plot_overlay(image, compare, title, thresh=None):
# shape and size of brain areas, we need to fix the alignment of the brains.
# The plot below shows that they are not yet aligned.

# load the freesurfer average brain
# load the subject's brain and the Freesurfer "fsaverage" template brain
subject_brain = nib.load(
op.join(misc_path, 'seeg', 'sample_seeg', 'mri', 'brain.mgz'))
template_brain = nib.load(
op.join(subjects_dir, 'fsaverage', 'mri', 'brain.mgz'))

Expand Down Expand Up @@ -290,13 +279,8 @@ def plot_overlay(image, compare, title, thresh=None):
# positions of all the voxels that had the contact's lookup number in
# the warped image.

# first we need our montage but it should be converted to "mri" coordinates
montage = raw.get_montage()
# ``subj_trans`` is "mri" to "head", we need "head" to "mri"
montage.apply_trans(mne.transforms.invert_transform(subj_trans))

montage_warped, elec_image, warped_elec_image = mne.warp_montage_volume(
montage, CT_aligned, reg_affine, sdr_morph,
elec_image, warped_elec_image = mne.warp_montage_volume(
raw.info, CT_aligned, reg_affine, sdr_morph,
subject_from='sample_seeg', subjects_dir=subjects_dir, thresh=CT_thresh)

fig, axes = plt.subplots(2, 1, figsize=(8, 8))
Expand All @@ -319,19 +303,10 @@ def plot_overlay(image, compare, title, thresh=None):

# sphinx_gallery_thumbnail_number = 8

# get native to head trans
fsaverage_trans = mne.channels.compute_native_head_t(montage_warped)

# set new montage
raw.set_montage(montage_warped)

# plot the resulting alignment
renderer = mne.viz.backends.renderer.create_3d_figure(**fig_kwargs)
brain = mne.viz.Brain('fsaverage', hemi='both', surf='pial',
cortex='low_contrast', alpha=0.2,
subjects_dir=subjects_dir, figure=renderer.figure)
brain = mne.viz.Brain('fsaverage', **brain_kwargs)
brain.add_sensors(raw)
mne.viz.set_3d_view(renderer.figure, **view_kwargs)
brain.show_view(**view_kwargs)

# %%
# This pipeline was developed based on previous work
Expand Down
14 changes: 2 additions & 12 deletions tutorials/io/30_reading_fnirs_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -245,19 +245,9 @@
subjects_dir = op.join(mne.datasets.sample.data_path(), 'subjects')
mne.datasets.fetch_fsaverage(subjects_dir=subjects_dir)

trans = mne.channels.compute_native_head_t(montage)

brain = mne.viz.Brain('fsaverage', subjects_dir=subjects_dir,
alpha=0.5, cortex='low_contrast')
# brain.add_head(dense=False)
brain.add_sensors(raw, trans)
brain.add_head(dense=False)
brain.add_sensors(raw)
brain.show_view(azimuth=90, elevation=90, distance=500,
focalpoint=(0., -10, 20))

fig = mne.viz.create_3d_figure(size=(800, 600), bgcolor='white')
fig = mne.viz.plot_alignment(
raw.info, trans=trans, subject='fsaverage', subjects_dir=subjects_dir,
surfaces=['brain', 'head'], coord_frame='mri', dig=True, show_axes=True,
fnirs=['channels', 'pairs', 'sources', 'detectors'], fig=fig)
mne.viz.set_3d_view(figure=fig, azimuth=90, elevation=90, distance=0.5,
focalpoint=(0., -0.01, 0.02))

0 comments on commit c1d78df

Please sign in to comment.