From b4c744b2a68531d212b4dd7eca5ec763d398ebf6 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 28 Jul 2021 12:53:04 -0700 Subject: [PATCH 1/3] refactor bem plotting code --- mne/_freesurfer.py | 68 +++++++++++++++++++++++++++++----------------- mne/report.py | 12 ++++---- mne/viz/misc.py | 55 ++++++++++++++----------------------- 3 files changed, 69 insertions(+), 66 deletions(-) diff --git a/mne/_freesurfer.py b/mne/_freesurfer.py index e0830d998d3..06e5fcb9564 100644 --- a/mne/_freesurfer.py +++ b/mne/_freesurfer.py @@ -57,44 +57,62 @@ def _import_nibabel(why='use MRI files'): return nib -def _mri_orientation(img, orientation): - """Get MRI orientation information from an image. - +def _reorient_image(img, axcodes='RAS'): + """Reorient an image to a given orientation. Parameters ---------- img : instance of SpatialImage The MRI image. + axcodes : tuple | str + The axis codes specifying the orientation, e.g. "RAS". + See :func:`nibabel.orientations.aff2axcodes`. + Returns + ------- + img_data : ndarray + The reoriented image data. + vox_ras_t : ndarray + The new transform from the new voxels to surface RAS. + Notes + ----- + .. versionadded:: 0.24 + """ + import nibabel as nib + orig_data = np.array(img.dataobj).astype(np.float32) + # reorient data to RAS + ornt = nib.orientations.axcodes2ornt( + nib.orientations.aff2axcodes(img.affine)).astype(int) + ras_ornt = nib.orientations.axcodes2ornt(axcodes) + ornt_trans = nib.orientations.ornt_transform(ornt, ras_ornt) + img_data = nib.orientations.apply_orientation(orig_data, ornt_trans) + orig_mgh = nib.MGHImage(orig_data, img.affine) + aff_trans = nib.orientations.inv_ornt_aff(ornt_trans, img.shape) + vox_ras_t = np.dot(orig_mgh.header.get_vox2ras_tkr(), aff_trans) + return img_data, vox_ras_t + + +def _mri_orientation(orientation): + """Get MRI orientation information from an image. + Parameters + ---------- orientation : str Orientation that you want. Can be "axial", "saggital", or "coronal". - Returns ------- - xyz : tuple, shape (3,) - The dimension indices for X, Y, and Z. - flips : tuple, shape (3,) - Whether each dimension requires a flip. - order : tuple, shape (3,) - The resulting order of the data if the given ``xyz`` and ``flips`` - are used. - + axis : int + The dimension of the axis to take slices over when plotting. + x : int + The dimension of the x axis. + y : int + The dimension of the y axis. Notes ----- .. versionadded:: 0.21 + .. versionchanged:: 0.24 """ - import nibabel as nib - _validate_type(img, nib.spatialimages.SpatialImage) _check_option('orientation', orientation, ('coronal', 'axial', 'sagittal')) - axcodes = ''.join(nib.orientations.aff2axcodes(img.affine)) - flips = {o: (1 if o in axcodes else -1) for o in 'RAS'} - axcodes = axcodes.replace('L', 'R').replace('P', 'A').replace('I', 'S') - order = dict( - coronal=('R', 'S', 'A'), - axial=('R', 'A', 'S'), - sagittal=('A', 'S', 'R'), - )[orientation] - xyz = tuple(axcodes.index(c) for c in order) - flips = tuple(flips[c] for c in order) - return xyz, flips, order + axis = dict(coronal=1, axial=2, sagittal=0)[orientation] + x, y = sorted(set([0, 1, 2]).difference(set([axis]))) + return axis, x, y def _get_mri_info_data(mri, data): diff --git a/mne/report.py b/mne/report.py index 6c79c9a2808..ce3e595b4ae 100644 --- a/mne/report.py +++ b/mne/report.py @@ -24,7 +24,7 @@ from .io import read_raw, read_info from .io._read_raw import supported as extension_reader_map from .io.pick import _DATA_CH_TYPES_SPLIT -from ._freesurfer import _mri_orientation +from ._freesurfer import _reorient_image, _mri_orientation from .utils import (logger, verbose, get_subjects_dir, warn, _ensure_int, fill_doc, _check_option, _validate_type, _safe_input) from .viz import (plot_events, plot_alignment, plot_cov, plot_projs_topomap, @@ -150,12 +150,9 @@ def _figs_to_mrislices(sl, n_jobs, **kwargs): parallel, p_fun, _ = parallel_func(_plot_mri_contours, use_jobs) outs = parallel(p_fun(slices=s, **kwargs) for s in np.array_split(sl, use_jobs)) - # deal with flip_z - flip_z = 1 out = list() for o in outs: - o, flip_z = o - out.extend(o[::flip_z]) + out.extend(o) return out @@ -1875,8 +1872,9 @@ def _render_one_bem_axis(self, mri_fname, surfaces, global_id, """Render one axis of bem contours (only PNG).""" import nibabel as nib nim = nib.load(mri_fname) - (_, _, z), _, _ = _mri_orientation(nim, orientation) - n_slices = nim.shape[z] + data = _reorient_image(nim)[0] + axis = _mri_orientation(orientation)[0] + n_slices = data.shape[axis] name = orientation html = [] diff --git a/mne/viz/misc.py b/mne/viz/misc.py index 824e261e89f..cdbfb3c0cdd 100644 --- a/mne/viz/misc.py +++ b/mne/viz/misc.py @@ -23,8 +23,8 @@ import numpy as np from ..defaults import DEFAULTS -from ..fixes import _get_img_fdata -from .._freesurfer import _mri_orientation, _read_mri_info, _check_mri +from .._freesurfer import (_reorient_image, _read_mri_info, _check_mri, + _mri_orientation) from ..rank import compute_rank from ..surface import read_surface from ..io.constants import FIFF @@ -32,7 +32,7 @@ from ..io.pick import (_DATA_CH_TYPES_SPLIT, pick_types, pick_info, pick_channels) from ..source_space import read_source_spaces, SourceSpaces, _ensure_src -from ..transforms import invert_transform, apply_trans, _frame_to_str +from ..transforms import apply_trans, _frame_to_str from ..utils import (logger, verbose, warn, _check_option, get_subjects_dir, _mask_to_onsets_offsets, _pl, _on_missing, fill_doc) from ..io.pick import _picks_by_type @@ -313,20 +313,14 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', _check_option('orientation', orientation, ('coronal', 'axial', 'sagittal')) # Load the T1 data - _, vox_mri_t, _, _, _, nim = _read_mri_info( + _, _, _, _, _, nim = _read_mri_info( mri_fname, units='mm', return_img=True) - mri_vox_t = invert_transform(vox_mri_t)['trans'] - del vox_mri_t - - # plot axes (x, y, z) as data axes - (x, y, z), (flip_x, flip_y, flip_z), order = _mri_orientation( - nim, orientation) - transpose = x < y - - data = _get_img_fdata(nim) - shift_x = data.shape[x] if flip_x < 0 else 0 - shift_y = data.shape[y] if flip_y < 0 else 0 - n_slices = data.shape[z] + + data, rasvox_mri_t = _reorient_image(nim) + mri_rasvox_t = np.linalg.inv(rasvox_mri_t) + axis, x, y = _mri_orientation(orientation) + + n_slices = data.shape[axis] if slices is None: slices = np.round(np.linspace(0, n_slices - 1, 14)).astype(int)[1:-1] slices = np.atleast_1d(slices).copy() @@ -337,9 +331,6 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', raise ValueError('slices must be a sorted 1D array of int with unique ' 'elements, at least one element, and no elements ' 'greater than %d, got %s' % (n_slices - 1, slices)) - if flip_z < 0: - # Proceed in the opposite order to maintain left-to-right / orientation - slices = slices[::-1] # create of list of surfaces surfs = list() @@ -347,7 +338,7 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', surf = dict() surf['rr'], surf['tris'] = read_surface(file_name) # move surface to voxel coordinate system - surf['rr'] = apply_trans(mri_vox_t, surf['rr']) + surf['rr'] = apply_trans(mri_rasvox_t, surf['rr']) surfs.append((surf, color)) sources = list() @@ -360,7 +351,7 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', f'{_frame_to_str[src[0]["coord_frame"]]}') for src_ in src: points = src_['rr'][src_['inuse'].astype(bool)] - sources.append(apply_trans(mri_vox_t, points * 1e3)) + sources.append(apply_trans(mri_rasvox_t, points * 1e3)) sources = np.concatenate(sources, axis=0) if img_output: @@ -384,17 +375,15 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', [[-np.inf], slices[:-1] + np.diff(slices) / 2., [np.inf]]) # float slicer = [slice(None)] * 3 ori_labels = dict(R='LR', A='PA', S='IS') - xlabels, ylabels = ori_labels[order[0]], ori_labels[order[1]] + xlabels, ylabels = ori_labels['RAS'[x]], ori_labels['RAS'[y]] path_effects = [patheffects.withStroke(linewidth=4, foreground="k", alpha=0.75)] out = list() if img_output else fig for ai, (ax, sl, lower, upper) in enumerate(zip( axs, slices, bounds[:-1], bounds[1:])): # adjust the orientations for good view - slicer[z] = sl - dat = data[tuple(slicer)] - dat = dat.T if transpose else dat - dat = dat[::flip_y, ::flip_x] + slicer[axis] = sl + dat = data[tuple(slicer)].T # First plot the anatomical data if img_output: @@ -408,16 +397,14 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', for surf, color in surfs: with warnings.catch_warnings(record=True): # ignore contour warn warnings.simplefilter('ignore') - ax.tricontour(flip_x * surf['rr'][:, x] + shift_x, - flip_y * surf['rr'][:, y] + shift_y, - surf['tris'], surf['rr'][:, z], + ax.tricontour(surf['rr'][:, x], surf['rr'][:, y], + surf['tris'], surf['rr'][:, axis], levels=[sl], colors=color, linewidths=1.0, zorder=1) if len(sources): - in_slice = (sources[:, z] >= lower) & (sources[:, z] < upper) - ax.scatter(flip_x * sources[in_slice, x] + shift_x, - flip_y * sources[in_slice, y] + shift_y, + in_slice = (sources[:, axis] >= lower) & (sources[:, axis] < upper) + ax.scatter(sources[in_slice, x], sources[in_slice, y], marker='.', color='#FF00FF', s=1, zorder=2) if show_indices: ax.text(dat.shape[1] // 8 + 0.5, 0.5, str(sl), @@ -448,7 +435,7 @@ def _plot_mri_contours(mri_fname, surfaces, src, orientation='coronal', fig.subplots_adjust(left=0., bottom=0., right=1., top=1., wspace=0., hspace=0.) plt_show(show, fig=fig) - return out, flip_z + return out def plot_bem(subject=None, subjects_dir=None, orientation='coronal', @@ -559,7 +546,7 @@ def plot_bem(subject=None, subjects_dir=None, orientation='coronal', # Plot the contours return _plot_mri_contours(mri_fname, surfaces, src, orientation, slices, - show, show_indices, show_orientation)[0] + show, show_indices, show_orientation) def _get_bem_plotting_surfaces(bem_path): From c7afecc430715201951e34aa0d32be05edc0c825 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 28 Jul 2021 13:14:12 -0700 Subject: [PATCH 2/3] fix docstyle --- mne/_freesurfer.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/mne/_freesurfer.py b/mne/_freesurfer.py index 06e5fcb9564..92b9cc3b6c0 100644 --- a/mne/_freesurfer.py +++ b/mne/_freesurfer.py @@ -59,6 +59,7 @@ def _import_nibabel(why='use MRI files'): def _reorient_image(img, axcodes='RAS'): """Reorient an image to a given orientation. + Parameters ---------- img : instance of SpatialImage @@ -66,12 +67,14 @@ def _reorient_image(img, axcodes='RAS'): axcodes : tuple | str The axis codes specifying the orientation, e.g. "RAS". See :func:`nibabel.orientations.aff2axcodes`. + Returns ------- img_data : ndarray The reoriented image data. vox_ras_t : ndarray The new transform from the new voxels to surface RAS. + Notes ----- .. versionadded:: 0.24 @@ -92,10 +95,12 @@ def _reorient_image(img, axcodes='RAS'): def _mri_orientation(orientation): """Get MRI orientation information from an image. + Parameters ---------- orientation : str Orientation that you want. Can be "axial", "saggital", or "coronal". + Returns ------- axis : int @@ -104,6 +109,7 @@ def _mri_orientation(orientation): The dimension of the x axis. y : int The dimension of the y axis. + Notes ----- .. versionadded:: 0.21 From 91214dbef78c3e745f7f7f15de620bd608d09deb Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 28 Jul 2021 13:47:04 -0700 Subject: [PATCH 3/3] clean up io in tutorials --- tutorials/forward/30_forward.py | 7 +++--- tutorials/intro/70_report.py | 41 +++++++++++++++++---------------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/tutorials/forward/30_forward.py b/tutorials/forward/30_forward.py index 1c4bb7399c4..bd8ee9678b8 100644 --- a/tutorials/forward/30_forward.py +++ b/tutorials/forward/30_forward.py @@ -18,9 +18,10 @@ data_path = sample.data_path() # the raw file containing the channel location + types -raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif' +sample_dir = op.join(data_path, 'MEG', 'sample',) +raw_fname = op.join(sample_dir, 'sample_audvis_raw.fif') # The paths to Freesurfer reconstructions -subjects_dir = data_path + '/subjects' +subjects_dir = op.join(data_path, 'subjects') subject = 'sample' # %% @@ -75,7 +76,7 @@ # alignment with the following code. # The transformation file obtained by coregistration -trans = data_path + '/MEG/sample/sample_audvis_raw-trans.fif' +trans = op.join(sample_dir, 'sample_audvis_raw-trans.fif') info = mne.io.read_info(raw_fname) # Here we look at the dense head, which isn't used for BEM computations but diff --git a/tutorials/intro/70_report.py b/tutorials/intro/70_report.py index 5d16a42d3d6..6340b6d2bc6 100644 --- a/tutorials/intro/70_report.py +++ b/tutorials/intro/70_report.py @@ -16,15 +16,19 @@ pages it generates are self-contained and do not require a running Python environment) but less flexible (you can't change code and re-run something directly within the browser). This tutorial covers the basics of building a -`~mne.Report`. As usual we'll start by importing the modules we need: +`~mne.Report`. As usual we'll start by importing the modules and data we need: """ # %% -import os +import os.path as op import matplotlib.pyplot as plt import mne +data_path = mne.datasets.sample.data_path(verbose=False) +sample_dir = op.join(data_path, 'MEG', 'sample') +subjects_dir = op.join(data_path, 'subjects') + # %% # Before getting started with :class:`mne.Report`, make sure the files you want # to render follow the filename conventions defined by MNE: @@ -77,9 +81,8 @@ # about not being able to render MRI and ``trans`` files without knowing the # subject. -path = mne.datasets.sample.data_path(verbose=False) report = mne.Report(verbose=True) -report.parse_folder(path, pattern='*raw.fif', render_bem=False) +report.parse_folder(data_path, pattern='*raw.fif', render_bem=False) report.save('report_basic.html', overwrite=True) # %% @@ -94,7 +97,7 @@ pattern = 'sample_audvis_filt-0-40_raw.fif' report = mne.Report(raw_psd=True, projs=True, verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_raw_psd.html', overwrite=True) # %% @@ -103,11 +106,10 @@ # containing an `~mne.Info` dictionary, from which the channel locations can be # read. -info_fname = os.path.join(path, 'MEG', 'sample', - 'sample_audvis_filt-0-40_raw.fif') +info_fname = op.join(sample_dir, 'sample_audvis_filt-0-40_raw.fif') pattern = 'sample_audvis_*proj.fif' report = mne.Report(info_fname=info_fname, verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_proj.html', overwrite=True) # %% @@ -118,9 +120,8 @@ # expensive, we'll also pass the ``mri_decim`` parameter for the benefit of our # documentation servers, and skip processing the :file:`.fif` files: -subjects_dir = os.path.join(path, 'subjects') report = mne.Report(subject='sample', subjects_dir=subjects_dir, verbose=True) -report.parse_folder(path, pattern='', mri_decim=25) +report.parse_folder(data_path, pattern='', mri_decim=25) report.save('report_mri_bem.html', overwrite=True) # %% @@ -131,7 +132,7 @@ pattern = 'sample_audvis-no-filter-ave.fif' report = mne.Report(verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_evoked.html', overwrite=True) # %% @@ -151,7 +152,7 @@ baseline = (None, 0) pattern = 'sample_audvis-no-filter-ave.fif' report = mne.Report(baseline=baseline, verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_evoked_baseline.html', overwrite=True) # %% @@ -160,10 +161,10 @@ # This will display ERP/ERF plots for both the original and whitened # :class:`~mne.Evoked` objects, but scalp topomaps only for the original. -cov_fname = os.path.join(path, 'MEG', 'sample', 'sample_audvis-cov.fif') +cov_fname = op.join(sample_dir, 'sample_audvis-cov.fif') baseline = (None, 0) report = mne.Report(cov_fname=cov_fname, baseline=baseline, verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_evoked_whitened.html', overwrite=True) # %% @@ -175,9 +176,9 @@ # information and should work): pattern = 'sample_audvis-cov.fif' -info_fname = os.path.join(path, 'MEG', 'sample', 'sample_audvis-ave.fif') +info_fname = op.join(sample_dir, 'sample_audvis-ave.fif') report = mne.Report(info_fname=info_fname, verbose=True) -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_cov.html', overwrite=True) # %% @@ -190,7 +191,7 @@ report = mne.Report(verbose=True) -fname_raw = os.path.join(path, 'MEG', 'sample', 'sample_audvis_raw.fif') +fname_raw = op.join(sample_dir, 'sample_audvis_raw.fif') raw = mne.io.read_raw_fif(fname_raw, verbose=False).crop(tmax=60) events = mne.find_events(raw, stim_channel='STI 014') event_id = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3, @@ -227,7 +228,7 @@ times = evoked_aud_left.times[::30] for t in times: figs.append(evoked_aud_left.plot_topomap(t, vmin=-300, vmax=300, res=100, - show=False)) + show=False)) plt.close(figs[-1]) report.add_slider_to_section(figs, times, 'Evoked Response', image_format='png') # can also use 'svg' @@ -245,7 +246,7 @@ report = mne.Report(info_fname=info_fname, subject='sample', subjects_dir=subjects_dir, verbose=True) pattern = "sample_audvis_raw-trans.fif" -report.parse_folder(path, pattern=pattern, render_bem=False) +report.parse_folder(data_path, pattern=pattern, render_bem=False) report.save('report_coreg.html', overwrite=True) # %% @@ -258,7 +259,7 @@ # :meth:`mne.viz.Brain.screenshot` method to save the figs in a slider. report = mne.Report(verbose=True) -fname_stc = os.path.join(path, 'MEG', 'sample', 'sample_audvis-meg') +fname_stc = op.join(sample_dir, 'sample_audvis-meg') stc = mne.read_source_estimate(fname_stc, subject='sample') figs = list() kwargs = dict(subjects_dir=subjects_dir, initial_time=0.13,