diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 4f04bfd8886..a4beb2caa01 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -15,6 +15,8 @@ Current (0.19.dev0) Changelog ~~~~~~~~~ +- Add :func:`mne.channels.make_dig_montage` to create :class:`mne.channels.DigMontage` objects out of np.arrays by `Joan Massich`_ + - Add support for reading in BrainVision CapTrak (BVCT) digitization coordinate files in :func:`mne.channels.read_dig_montage` by `Stefan Appelhoff`_ - Allow :meth:`mne.Annotations.crop` to support negative ``tmin`` and ``tmax`` by `Joan Massich`_ diff --git a/doc/conf.py b/doc/conf.py index 634e66faaff..bdceb5c5739 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -540,6 +540,7 @@ def reset_warnings(gallery_conf, fname): 'Covariance': 'mne.Covariance', 'Annotations': 'mne.Annotations', 'Montage': 'mne.channels.Montage', 'DigMontage': 'mne.channels.DigMontage', + 'Digitization': 'mne.digitization.Digitization', 'VectorSourceEstimate': 'mne.VectorSourceEstimate', 'VolSourceEstimate': 'mne.VolSourceEstimate', 'VolVectorSourceEstimate': 'mne.VolVectorSourceEstimate', diff --git a/doc/python_reference.rst b/doc/python_reference.rst index 6796486e806..5ced5e9c709 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -310,6 +310,7 @@ Projections: read_montage get_builtin_montages read_dig_montage + make_dig_montage read_layout find_layout make_eeg_layout diff --git a/examples/visualization/plot_3d_to_2d.py b/examples/visualization/plot_3d_to_2d.py index 43ef9a22771..ec595e04420 100644 --- a/examples/visualization/plot_3d_to_2d.py +++ b/examples/visualization/plot_3d_to_2d.py @@ -50,9 +50,8 @@ mat = loadmat(path_data) ch_names = mat['ch_names'].tolist() elec = mat['elec'] # electrode coordinates in meters -dig_ch_pos = dict(zip(ch_names, elec)) -mon = mne.channels.DigMontage(dig_ch_pos=dig_ch_pos) -info = mne.create_info(ch_names, 1000., 'ecog', montage=mon) +montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec))) +info = mne.create_info(ch_names, 1000., 'ecog', montage=montage) print('Created %s channel positions' % len(ch_names)) ############################################################################### @@ -67,7 +66,7 @@ fig = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, surfaces=['pial'], meg=False) set_3d_view(figure=fig, azimuth=200, elevation=70) -xy, im = snapshot_brain_montage(fig, mon) +xy, im = snapshot_brain_montage(fig, montage) # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in info['ch_names']]) diff --git a/mne/channels/__init__.py b/mne/channels/__init__.py index f0bb522cbb6..dda46f0c990 100644 --- a/mne/channels/__init__.py +++ b/mne/channels/__init__.py @@ -6,7 +6,7 @@ from .layout import (Layout, make_eeg_layout, make_grid_layout, read_layout, find_layout, generate_2d_layout) from .montage import (read_montage, read_dig_montage, Montage, DigMontage, - get_builtin_montages) + get_builtin_montages, make_dig_montage) from .channels import (equalize_channels, rename_channels, fix_mag_coil_types, read_ch_connectivity, _get_ch_type, find_ch_connectivity, make_1020_channel_selections) diff --git a/mne/channels/_dig_montage_utils.py b/mne/channels/_dig_montage_utils.py new file mode 100644 index 00000000000..bb1dbd9edde --- /dev/null +++ b/mne/channels/_dig_montage_utils.py @@ -0,0 +1,306 @@ +# Authors: Alexandre Gramfort +# Denis Engemann +# Martin Luessi +# Eric Larson +# Marijn van Vliet +# Jona Sassenhagen +# Teon Brooks +# Christian Brodbeck +# Stefan Appelhoff +# Joan Massich +# +# License: Simplified BSD + +import xml.etree.ElementTree as ElementTree + +import numpy as np + +from ..transforms import apply_trans, get_ras_to_neuromag_trans + +from ..io.constants import FIFF +from ..io.open import fiff_open +from ..digitization._utils import _read_dig_fif +from ..utils import _check_fname, Bunch, warn + + +def _fix_data_fiducials(data): + nasion, rpa, lpa = data.nasion, data.rpa, data.lpa + if any(x is None for x in (nasion, rpa, lpa)): + if data.elp is None or data.point_names is None: + raise ValueError('ELP points and names must be specified for ' + 'transformation.') + names = [name.lower() for name in data.point_names] + + # check that all needed points are present + kinds = ('nasion', 'lpa', 'rpa') + missing = [name for name in kinds if name not in names] + if len(missing) > 0: + raise ValueError('The points %s are missing, but are needed ' + 'to transform the points to the MNE ' + 'coordinate system. Either add the points, ' + 'or read the montage with transform=False.' + % str(missing)) + + data.nasion, data.lpa, data.rpa = [ + data.elp[names.index(kind)] for kind in kinds + ] + + # remove fiducials from elp + mask = np.ones(len(names), dtype=bool) + for fid in ['nasion', 'lpa', 'rpa']: + mask[names.index(fid)] = False + data.elp = data.elp[mask] + data.point_names = [p for pi, p in enumerate(data.point_names) + if mask[pi]] + return data + + +def _transform_to_head_call(data): + """Transform digitizer points to Neuromag head coordinates. + + Parameters + ---------- + data : Bunch. + replicates DigMontage old structure. Requires the following fields: + ['nasion', 'lpa', 'rpa', 'hsp', 'hpi', 'elp', 'coord_frame', + 'dig_ch_pos'] + + Returns + ------- + data : Bunch. + transformed version of input data. + """ + if data.coord_frame == 'head': # nothing to do + return data + nasion, rpa, lpa = data.nasion, data.rpa, data.lpa + + native_head_t = get_ras_to_neuromag_trans(nasion, lpa, rpa) + data.nasion, data.lpa, data.rpa = apply_trans( + native_head_t, np.array([nasion, lpa, rpa])) + if data.elp is not None: + data.elp = apply_trans(native_head_t, data.elp) + if data.hsp is not None: + data.hsp = apply_trans(native_head_t, data.hsp) + if data.dig_ch_pos is not None: + for key, val in data.dig_ch_pos.items(): + data.dig_ch_pos[key] = apply_trans(native_head_t, val) + data.coord_frame = 'head' + + return data + + +_cardinal_ident_mapping = { + FIFF.FIFFV_POINT_NASION: 'nasion', + FIFF.FIFFV_POINT_LPA: 'lpa', + FIFF.FIFFV_POINT_RPA: 'rpa', +} + + +def _read_dig_montage_fif( + fname, + _raise_transform_err, + _all_data_kwargs_are_none, +): + from .montage import _check_frame # circular dep + + if _raise_transform_err: + raise ValueError('transform must be True and dev_head_t must be ' + 'False for FIF dig montage') + if not _all_data_kwargs_are_none: + raise ValueError('hsp, hpi, elp, point_names, egi must all be ' + 'None if fif is not None') + + _check_fname(fname, overwrite='read', must_exist=True) + # Load the dig data + f, tree = fiff_open(fname)[:2] + with f as fid: + dig = _read_dig_fif(fid, tree) + + # Split up the dig points by category + hsp = list() + hpi = list() + elp = list() + point_names = list() + fids = dict() + dig_ch_pos = dict() + for d in dig: + if d['kind'] == FIFF.FIFFV_POINT_CARDINAL: + _check_frame(d, 'head') + fids[_cardinal_ident_mapping[d['ident']]] = d['r'] + elif d['kind'] == FIFF.FIFFV_POINT_HPI: + _check_frame(d, 'head') + hpi.append(d['r']) + elp.append(d['r']) + point_names.append('HPI%03d' % d['ident']) + elif d['kind'] == FIFF.FIFFV_POINT_EXTRA: + _check_frame(d, 'head') + hsp.append(d['r']) + elif d['kind'] == FIFF.FIFFV_POINT_EEG: + _check_frame(d, 'head') + dig_ch_pos['EEG%03d' % d['ident']] = d['r'] + + return Bunch( + nasion=fids['nasion'], lpa=fids['lpa'], rpa=fids['rpa'], + hsp=np.array(hsp) if len(hsp) else None, + hpi=hpi, # XXX: why this is returned as empty list instead of None? + elp=np.array(elp) if len(elp) else None, + point_names=point_names, + dig_ch_pos=dig_ch_pos, + coord_frame='head', + ) + + +def _read_dig_montage_egi( + fname, + _scaling, + _all_data_kwargs_are_none, +): + + if not _all_data_kwargs_are_none: + raise ValueError('hsp, hpi, elp, point_names, fif must all be ' + 'None if egi is not None') + _check_fname(fname, overwrite='read', must_exist=True) + + root = ElementTree.parse(fname).getroot() + ns = root.tag[root.tag.index('{'):root.tag.index('}') + 1] + sensors = root.find('%ssensorLayout/%ssensors' % (ns, ns)) + fids = dict() + dig_ch_pos = dict() + + fid_name_map = {'Nasion': 'nasion', + 'Right periauricular point': 'rpa', + 'Left periauricular point': 'lpa'} + + for s in sensors: + name, number, kind = s[0].text, int(s[1].text), int(s[2].text) + coordinates = np.array([float(s[3].text), float(s[4].text), + float(s[5].text)]) + + coordinates *= _scaling + + # EEG Channels + if kind == 0: + dig_ch_pos['EEG %03d' % number] = coordinates + # Reference + elif kind == 1: + dig_ch_pos['EEG %03d' % + (len(dig_ch_pos.keys()) + 1)] = coordinates + # XXX: we should do something with this (ref and eeg get mixed) + + # Fiducials + elif kind == 2: + fid_name = fid_name_map[name] + fids[fid_name] = coordinates + # Unknown + else: + warn('Unknown sensor type %s detected. Skipping sensor...' + 'Proceed with caution!' % kind) + + return Bunch( + # EGI stuff + nasion=fids['nasion'], lpa=fids['lpa'], rpa=fids['rpa'], + dig_ch_pos=dig_ch_pos, coord_frame='unknown', + # not EGI stuff + hsp=None, hpi=None, elp=None, point_names=None, + ) + + +def _foo_get_data_from_dig(dig): + # XXXX: + # This does something really similar to _read_dig_montage_fif but: + # - does not check coord_frame + # - does not do any operation that implies assumptions with the names + + # Split up the dig points by category + hsp, hpi, elp = list(), list(), list() + fids, dig_ch_pos_location = dict(), list() + + for d in dig: + if d['kind'] == FIFF.FIFFV_POINT_CARDINAL: + fids[_cardinal_ident_mapping[d['ident']]] = d['r'] + elif d['kind'] == FIFF.FIFFV_POINT_HPI: + hpi.append(d['r']) + elp.append(d['r']) + # XXX: point_names.append('HPI%03d' % d['ident']) + elif d['kind'] == FIFF.FIFFV_POINT_EXTRA: + hsp.append(d['r']) + elif d['kind'] == FIFF.FIFFV_POINT_EEG: + # XXX: dig_ch_pos['EEG%03d' % d['ident']] = d['r'] + dig_ch_pos_location.append(d['r']) + + dig_coord_frames = set([d['coord_frame'] for d in dig]) + assert len(dig_coord_frames) == 1, 'Only single coordinate frame in dig is supported' # noqa # XXX + + return Bunch( + nasion=fids.get('nasion', None), + lpa=fids.get('lpa', None), + rpa=fids.get('rpa', None), + hsp=np.array(hsp) if len(hsp) else None, + hpi=np.array(hpi) if len(hpi) else None, + elp=np.array(elp) if len(elp) else None, + dig_ch_pos_location=dig_ch_pos_location, + coord_frame=dig_coord_frames.pop(), + ) + + +def _read_dig_montage_bvct( + fname, + unit, + _all_data_kwargs_are_none, +): + + if not _all_data_kwargs_are_none: + raise ValueError('hsp, hpi, elp, point_names, fif must all be ' + 'None if egi is not None') + _check_fname(fname, overwrite='read', must_exist=True) + + # CapTrak is natively in mm + scale = dict(mm=1e-3, cm=1e-2, auto=1e-3, m=1) + if unit not in scale: + raise ValueError("Unit needs to be one of %s, not %r" % + (sorted(scale.keys()), unit)) + if unit not in ['mm', 'auto']: + warn('Using "{}" as unit for BVCT file. BVCT files are usually ' + 'specified in "mm". This might lead to errors.'.format(unit), + RuntimeWarning) + + root = ElementTree.parse(fname).getroot() + sensors = root.find('CapTrakElectrodeList') + + fids = {} + dig_ch_pos = {} + + fid_name_map = {'Nasion': 'nasion', 'RPA': 'rpa', 'LPA': 'lpa'} + + for s in sensors: + name = s.find('Name').text + + # Need to prune "GND" and "REF": these are not included in the raw + # data and will raise errors when we try to do raw.set_montage(...) + # XXX eventually this should be stored in ch['loc'][3:6] + # but we don't currently have such capabilities here + if name in ['GND', 'REF']: + continue + + fid = name in fid_name_map + coordinates = np.array([float(s.find('X').text), + float(s.find('Y').text), + float(s.find('Z').text)]) + + coordinates *= scale[unit] + + # Fiducials + if fid: + fid_name = fid_name_map[name] + fids[fid_name] = coordinates + # EEG Channels + else: + dig_ch_pos[name] = coordinates + + return Bunch( + # BVCT stuff + nasion=fids['nasion'], lpa=fids['lpa'], rpa=fids['rpa'], + dig_ch_pos=dig_ch_pos, coord_frame='unknown', + # not BVCT stuff + hsp=None, hpi=None, elp=None, point_names=None, + ) diff --git a/mne/channels/montage.py b/mne/channels/montage.py index 43d2689ed25..f925bd6c986 100644 --- a/mne/channels/montage.py +++ b/mne/channels/montage.py @@ -7,6 +7,7 @@ # Teon Brooks # Christian Brodbeck # Stefan Appelhoff +# Joan Massich # # License: Simplified BSD @@ -21,15 +22,28 @@ from .channels import _contains_ch_type from ..transforms import (apply_trans, get_ras_to_neuromag_trans, _sph_to_cart, _topo_to_sph, _str_to_frame, _frame_to_str) +from ..digitization import Digitization from ..digitization._utils import (_make_dig_points, _read_dig_points, - _read_dig_fif, write_dig) + write_dig) from ..io.pick import pick_types -from ..io.open import fiff_open from ..io.constants import FIFF -from ..utils import (_check_fname, warn, copy_function_doc_to_method_doc, - _check_option) +from ..utils import (warn, copy_function_doc_to_method_doc, + _check_option, Bunch, deprecated, _validate_type) from .layout import _pol_to_cart, _cart_to_sph +from ._dig_montage_utils import _transform_to_head_call, _read_dig_montage_fif +from ._dig_montage_utils import _read_dig_montage_egi, _read_dig_montage_bvct +from ._dig_montage_utils import _foo_get_data_from_dig +from ._dig_montage_utils import _fix_data_fiducials + +DEPRECATED_PARAM = object() + + +def _check_get_coord_frame(dig): + _MSG = 'Only single coordinate frame in dig is supported' + dig_coord_frames = set([d['coord_frame'] for d in dig]) + assert len(dig_coord_frames) == 1, _MSG + return _frame_to_str[dig_coord_frames.pop()] class Montage(object): @@ -394,6 +408,98 @@ def read_montage(kind, ch_names=None, path=None, unit='m', transform=False): lpa=fids['lpa'], nasion=fids['nasion'], rpa=fids['rpa']) +def make_dig_montage(ch_pos=None, nasion=None, lpa=None, rpa=None, + hsp=None, hpi=None, hpi_dev=None, coord_frame='unknown', + transform_to_head=False, compute_dev_head_t=False): + r"""Make montage from arrays. + + Parameters + ---------- + ch_pos : dict + Dictionary of channel positions. Keys are channel names and values + are 3D coordinates - array of shape (3,) - in native digitizer space + in m. + nasion : None | array, shape (3,) + The position of the nasion fiducial point. + This point is assumed to be in the native digitizer space in m. + lpa : None | array, shape (3,) + The position of the left periauricular fiducial point. + This point is assumed to be in the native digitizer space in m. + rpa : None | array, shape (3,) + The position of the right periauricular fiducial point. + This point is assumed to be in the native digitizer space in m. + hsp : None | array, shape (n_points, 3) + This corresponds to an array of positions of the headshape points in + 3d. These points are assumed to be in the native digitizer space in m. + hpi : None | array, shape (n_hpi, 3) + This corresponds to an array of HPI points in the native digitizer + space. They only necessary if computation of a ``compute_dev_head_t`` + is True. + hpi_dev : None | array, shape (n_hpi, 3) + This corresponds to an array of HPI points. These points are in device + space, and are only necessary if computation of a + ``compute_dev_head_t`` is True. + coord_frame : str + The coordinate frame of the points. Usually this is "unknown" + for native digitizer space. + transform_to_head : bool + If True (default), points will be transformed to Neuromag head space. + The fiducials (nasion, lpa, and rpa) must be specified. This is useful + for points captured using a device that does not automatically convert + points to Neuromag head coordinates + (e.g., Polhemus FastSCAN). + compute_dev_head_t : bool + If True, a Dev-to-Head transformation matrix will be added to the + montage. To get a proper `dev_head_t`, the hpi and the hpi_dev points + must be in the same order. If False (default), no transformation will + be added to the montage. + + Returns + ------- + montage : instance of DigMontage + The montage. + + See Also + -------- + Montage + read_montage + DigMontage + read_dig_montage + + """ + # XXX: hpi was historically elp + # XXX: hpi_dev was historically hpi + assert coord_frame in ('unknown', 'head') + from ..coreg import fit_matched_points + data = Bunch( + nasion=nasion, lpa=lpa, rpa=rpa, + elp=hpi, dig_ch_pos=ch_pos, hsp=hsp, + coord_frame=coord_frame, + ) + if transform_to_head: + data = _transform_to_head_call(data) + + if compute_dev_head_t: + if data.elp is None or hpi_dev is None: + raise RuntimeError('must have both elp and hpi to compute the ' + 'device to head transform') + else: + # here is hpi + dev_head_t = fit_matched_points( + tgt_pts=data.elp, src_pts=hpi_dev, out='trans' + ) # XXX: shall we make it a Transform? rather than np.array + else: + dev_head_t = None + + ch_names = list() if ch_pos is None else list(sorted(ch_pos.keys())) + dig = _make_dig_points( + nasion=data.nasion, lpa=data.lpa, rpa=data.rpa, hpi=data.elp, + extra_points=data.hsp, dig_ch_pos=data.dig_ch_pos, + coord_frame=data.coord_frame, + ) + return DigMontage(dig=dig, ch_names=ch_names, dev_head_t=dev_head_t) + + class DigMontage(object): """Montage for digitized electrode and headshape position data. @@ -406,24 +512,32 @@ class DigMontage(object): hsp : array, shape (n_points, 3) The positions of the headshape points in 3d. These points are in the native digitizer space. + Deprecated, will be removed in 0.20. hpi : array, shape (n_hpi, 3) The positions of the head-position indicator coils in 3d. These points are in the MEG device space. + Deprecated, will be removed in 0.20. elp : array, shape (n_hpi, 3) The positions of the head-position indicator coils in 3d. This is typically in the native digitizer space. + Deprecated, will be removed in 0.20. point_names : list, shape (n_elp) The names of the digitized points for hpi and elp. + Deprecated, will be removed in 0.20. nasion : array, shape (1, 3) - The position of the nasion fidicual point. + The position of the nasion fiducial point. + Deprecated, will be removed in 0.20. lpa : array, shape (1, 3) - The position of the left periauricular fidicual point. + The position of the left periauricular fiducial point. + Deprecated, will be removed in 0.20. rpa : array, shape (1, 3) - The position of the right periauricular fidicual point. + The position of the right periauricular fiducial point. + Deprecated, will be removed in 0.20. dev_head_t : array, shape (4, 4) A Device-to-Head transformation matrix. dig_ch_pos : dict Dictionary of channel positions given in meters. + Deprecated, will be removed in 0.20. .. versionadded:: 0.12 @@ -431,6 +545,13 @@ class DigMontage(object): The coordinate frame of the points. Usually this is "unknown" for native digitizer space. + .. versionadded:: 0.19 + + dig : instance of Digitization + The object containing all the dig points. + ch_names : list of str + The names of the EEG channels. + See Also -------- Montage @@ -442,107 +563,140 @@ class DigMontage(object): .. versionadded:: 0.9.0 """ - def __init__(self, hsp=None, hpi=None, elp=None, point_names=None, - nasion=None, lpa=None, rpa=None, dev_head_t=None, - dig_ch_pos=None, coord_frame='unknown'): # noqa: D102 - self.hsp = hsp - self.hpi = hpi - if elp is not None: - if not isinstance(point_names, Iterable): - raise TypeError('If elp is specified, point_names must ' - 'provide a list of str with one entry per ELP ' - 'point') - point_names = list(point_names) - if len(point_names) != len(elp): - raise ValueError('elp contains %i points but %i ' - 'point_names were specified.' % - (len(elp), len(point_names))) - self.elp = elp - self.point_names = point_names - - self.nasion = nasion - self.lpa = lpa - self.rpa = rpa - self.dev_head_t = dev_head_t - self.dig_ch_pos = dig_ch_pos - if not isinstance(coord_frame, str) or \ - coord_frame not in _str_to_frame: - raise ValueError('coord_frame must be one of %s, got %s' - % (sorted(_str_to_frame.keys()), coord_frame)) - self.coord_frame = coord_frame + def __init__(self, + hsp=DEPRECATED_PARAM, hpi=DEPRECATED_PARAM, elp=DEPRECATED_PARAM, + point_names=DEPRECATED_PARAM, nasion=DEPRECATED_PARAM, + lpa=DEPRECATED_PARAM, rpa=DEPRECATED_PARAM, + dev_head_t=None, dig_ch_pos=DEPRECATED_PARAM, + coord_frame=DEPRECATED_PARAM, + dig=None, ch_names=None, + ): # noqa: D102 + # XXX: dev_head_t now is np.array, we should add dev_head_transform + # (being instance of Transformation) and move the parameter to the + # end of the call. + _non_deprecated_kwargs = [ + key for key, val in dict( + hsp=hsp, hpi=hpi, elp=elp, point_names=point_names, + nasion=nasion, lpa=lpa, rpa=rpa, + dig_ch_pos=dig_ch_pos, coord_frame=coord_frame, + ).items() if val is not DEPRECATED_PARAM + ] + if not _non_deprecated_kwargs: + _validate_type(item=dig, types=Digitization, + item_name='dig', type_name='Digitization') + ch_names = list() if ch_names is None else ch_names + n_eeg = sum([1 for d in dig if d['kind'] == FIFF.FIFFV_POINT_EEG]) + if n_eeg != len(ch_names): + raise ValueError( + 'The number of EEG channels (%d) does not match the number' + ' of channel names provided (%d)' % (n_eeg, len(ch_names)) + ) + + self.dev_head_t = dev_head_t + self.dig = dig + self.ch_names = ch_names + self._coord_frame = _check_get_coord_frame(self.dig) + else: + # Deprecated + _msg = ( + "Using {params} in DigMontage constructor is deprecated." + " Use 'dig', and 'ch_names' instead." + ).format(params=", ".join( + ["'{}'".format(k) for k in _non_deprecated_kwargs] + )) + warn(_msg, DeprecationWarning) + + # Restore old defaults + hsp = None if hsp is DEPRECATED_PARAM else hsp + hpi = None if hpi is DEPRECATED_PARAM else hpi + elp = None if elp is DEPRECATED_PARAM else elp + nasion = None if nasion is DEPRECATED_PARAM else nasion + lpa = None if lpa is DEPRECATED_PARAM else lpa + rpa = None if rpa is DEPRECATED_PARAM else rpa + dig_ch_pos = None if dig_ch_pos is DEPRECATED_PARAM else dig_ch_pos + coord_frame = \ + 'unknown' if coord_frame is DEPRECATED_PARAM else coord_frame + point_names = \ + None if point_names is DEPRECATED_PARAM else point_names + + # Old behavior + if elp is not None: + if not isinstance(point_names, Iterable): + raise TypeError('If elp is specified, point_names must' + ' provide a list of str with one entry per' + ' ELP point.') + point_names = list(point_names) + if len(point_names) != len(elp): + raise ValueError('elp contains %i points but %i ' + 'point_names were specified.' % + (len(elp), len(point_names))) + + self.dev_head_t = dev_head_t + self._point_names = point_names + self.ch_names = \ + [] if dig_ch_pos is None else list(sorted(dig_ch_pos.keys())) + self._hpi = hpi + self._coord_frame = coord_frame + self.dig = _make_dig_points( + nasion=nasion, lpa=lpa, rpa=rpa, hpi=elp, + extra_points=hsp, dig_ch_pos=dig_ch_pos, + coord_frame=self._coord_frame, + ) + + @property + def point_names(self): + warn('"point_names" attribute is deprecated and will be removed' + ' in v0.20', DeprecationWarning) + return self._point_names + + @property + def coord_frame(self): + warn('"coord_frame" attribute is deprecated and will be removed' + ' in v0.20', DeprecationWarning) + return self._coord_frame def __repr__(self): """Return string representation.""" + _data = _foo_get_data_from_dig(self.dig) s = ('' % - (len(self.hsp) if self.hsp is not None else 0, - len(self.point_names) if self.point_names is not None else 0, - sum(x is not None for x in (self.lpa, self.rpa, self.nasion)), - len(self.dig_ch_pos) if self.dig_ch_pos is not None else 0,)) + (len(_data.hsp) if _data.hsp is not None else 0, + len(_data.hpi) if _data.hpi is not None else 0, + sum(x is not None for x in (_data.lpa, _data.rpa, _data.nasion)), + len(_data.dig_ch_pos_location) if _data.dig_ch_pos_location is not None else 0,)) # noqa return s @copy_function_doc_to_method_doc(plot_montage) def plot(self, scale_factor=20, show_names=False, kind='3d', show=True): + # XXX: plot_montage takes an empty info and sets 'self' + # Therefore it should not be a representation problem. return plot_montage(self, scale_factor=scale_factor, show_names=show_names, kind=kind, show=show) def transform_to_head(self): """Transform digitizer points to Neuromag head coordinates.""" - if self.coord_frame == 'head': # nothing to do - return - nasion, rpa, lpa = self.nasion, self.rpa, self.lpa - if any(x is None for x in (nasion, rpa, lpa)): - if self.elp is None or self.point_names is None: - raise ValueError('ELP points and names must be specified for ' - 'transformation.') - names = [name.lower() for name in self.point_names] - - # check that all needed points are present - kinds = ('nasion', 'lpa', 'rpa') - missing = [name for name in kinds if name not in names] - if len(missing) > 0: - raise ValueError('The points %s are missing, but are needed ' - 'to transform the points to the MNE ' - 'coordinate system. Either add the points, ' - 'or read the montage with transform=False.' - % str(missing)) - - nasion, lpa, rpa = [self.elp[names.index(kind)] for kind in kinds] - - # remove fiducials from elp - mask = np.ones(len(names), dtype=bool) - for fid in ['nasion', 'lpa', 'rpa']: - mask[names.index(fid)] = False - self.elp = self.elp[mask] - self.point_names = [p for pi, p in enumerate(self.point_names) - if mask[pi]] - - native_head_t = get_ras_to_neuromag_trans(nasion, lpa, rpa) - self.nasion, self.lpa, self.rpa = apply_trans( - native_head_t, np.array([nasion, lpa, rpa])) - if self.elp is not None: - self.elp = apply_trans(native_head_t, self.elp) - if self.hsp is not None: - self.hsp = apply_trans(native_head_t, self.hsp) - if self.dig_ch_pos is not None: - for key, val in self.dig_ch_pos.items(): - self.dig_ch_pos[key] = apply_trans(native_head_t, val) - self.coord_frame = 'head' - + raise RuntimeError('The transform_to_head method has been removed to ' + 'enforce that DigMontage are constructed already ' + 'in the correct coordinate system. This method ' + 'will disappear in version 0.20.') + + @deprecated( + 'compute_dev_head_t is deprecated and will be removed in 0.20.' + ) def compute_dev_head_t(self): """Compute the Neuromag dev_head_t from matched points.""" + if not hasattr(self, '_hpi'): + raise RuntimeError( + 'Cannot compute dev_head_t if DigMontage was not created' + ' from arrays') + from ..coreg import fit_matched_points - if self.elp is None or self.hpi is None: + data = _foo_get_data_from_dig(self.dig) + if data.elp is None or self._hpi is None: raise RuntimeError('must have both elp and hpi to compute the ' 'device to head transform') - self.dev_head_t = fit_matched_points(tgt_pts=self.elp, - src_pts=self.hpi, out='trans') - - def _get_dig(self): - """Get the digitization list.""" - return _make_dig_points( - nasion=self.nasion, lpa=self.lpa, rpa=self.rpa, hpi=self.elp, - extra_points=self.hsp, dig_ch_pos=self.dig_ch_pos) + self.dev_head_t = fit_matched_points(tgt_pts=data.elp, + src_pts=self._hpi, out='trans') def save(self, fname): """Save digitization points to FIF. @@ -552,17 +706,56 @@ def save(self, fname): fname : str The filename to use. Should end in .fif or .fif.gz. """ - if self.coord_frame != 'head': + if self._coord_frame != 'head': raise RuntimeError('Can only write out digitization points in ' 'head coordinates.') - write_dig(fname, self._get_dig()) - - -_cardinal_ident_mapping = { - FIFF.FIFFV_POINT_NASION: 'nasion', - FIFF.FIFFV_POINT_LPA: 'lpa', - FIFF.FIFFV_POINT_RPA: 'rpa', -} + write_dig(fname, self.dig) + + @property + def dig_ch_pos(self): + warn('"dig_ch_pos" attribute is deprecated and will be removed in ' + 'v0.20', DeprecationWarning) + return self._ch_pos() + + def _get_ch_pos(self): + return dict(zip(self.ch_names, + _foo_get_data_from_dig(self.dig).dig_ch_pos_location)) + + @property + def elp(self): + warn('"elp" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return _foo_get_data_from_dig(self.dig).elp + + @property + def hpi(self): + warn('"hpi" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return getattr(self, '_hpi', None) + + @property + def hsp(self): + warn('"hsp" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return _foo_get_data_from_dig(self.dig).hsp + + @property + def lpa(self): + warn('"lpa" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return _foo_get_data_from_dig(self.dig).lpa + + @property + def rpa(self): + warn('"rpa" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return _foo_get_data_from_dig(self.dig).rpa + + @property + def nasion(self): + warn('"nasion" attribute is deprecated and will be removed in v0.20', + DeprecationWarning) + return _foo_get_data_from_dig(self.dig).nasion def _check_frame(d, frame_str): @@ -572,9 +765,17 @@ def _check_frame(d, frame_str): % (frame_str, _frame_to_str[d['coord_frame']])) -def read_dig_montage(hsp=None, hpi=None, elp=None, point_names=None, - unit='auto', fif=None, egi=None, bvct=None, - transform=True, dev_head_t=False, ): +def _get_scaling(unit, scale): + if unit not in scale: + raise ValueError("Unit needs to be one of %s, not %r" % + (sorted(scale.keys()), unit)) + else: + return scale[unit] + + +def read_dig_montage(hsp=None, hpi=None, elp=None, + point_names=None, unit='auto', fif=None, egi=None, + bvct=None, transform=True, dev_head_t=False, ): r"""Read subject-specific digitization montage from a file. Parameters @@ -629,7 +830,7 @@ def read_dig_montage(hsp=None, hpi=None, elp=None, point_names=None, transform : bool If True (default), points will be transformed to Neuromag space using :meth:`DigMontage.transform_to_head`. - The fidicuals (nasion, lpa, and rpa) must be specified. + The fiducials (nasion, lpa, and rpa) must be specified. This is useful for points captured using a device that does not automatically convert points to Neuromag head coordinates (e.g., Polhemus FastSCAN). @@ -658,157 +859,43 @@ def read_dig_montage(hsp=None, hpi=None, elp=None, point_names=None, .. versionadded:: 0.9.0 """ + # XXX: This scaling business seems really dangerous to me. + EGI_SCALE = dict(mm=1e-3, cm=1e-2, auto=1e-2, m=1) + NUMPY_DATA_SCALE = dict(mm=1e-3, cm=1e-2, auto=1e-3, m=1) + if fif is not None: - # Use a different code path - if dev_head_t or not transform: - raise ValueError('transform must be True and dev_head_t must be ' - 'False for FIF dig montage') - if not all(x is None for x in (hsp, hpi, elp, point_names, egi, bvct)): - raise ValueError('hsp, hpi, elp, point_names, egi, bvct must all ' - 'be None if fif is not None.') - _check_fname(fif, overwrite='read', must_exist=True) - # Load the dig data - f, tree = fiff_open(fif)[:2] - with f as fid: - dig = _read_dig_fif(fid, tree) - # Split up the dig points by category - hsp = list() - hpi = list() - elp = list() - point_names = list() - fids = dict() - dig_ch_pos = dict() - for d in dig: - if d['kind'] == FIFF.FIFFV_POINT_CARDINAL: - _check_frame(d, 'head') - fids[_cardinal_ident_mapping[d['ident']]] = d['r'] - elif d['kind'] == FIFF.FIFFV_POINT_HPI: - _check_frame(d, 'head') - hpi.append(d['r']) - elp.append(d['r']) - point_names.append('HPI%03d' % d['ident']) - elif d['kind'] == FIFF.FIFFV_POINT_EXTRA: - _check_frame(d, 'head') - hsp.append(d['r']) - elif d['kind'] == FIFF.FIFFV_POINT_EEG: - _check_frame(d, 'head') - dig_ch_pos['EEG%03d' % d['ident']] = d['r'] - fids = [fids.get(key) for key in ('nasion', 'lpa', 'rpa')] - hsp = np.array(hsp) if len(hsp) else None - elp = np.array(elp) if len(elp) else None - coord_frame = 'head' + _raise_transform_err = True if dev_head_t or not transform else False + data = _read_dig_montage_fif( + fname=fif, + _raise_transform_err=_raise_transform_err, + _all_data_kwargs_are_none=all( + x is None for x in (hsp, hpi, elp, point_names, egi, bvct)) + ) elif egi is not None: - if not all(x is None for x in (hsp, hpi, elp, point_names, fif, bvct)): - raise ValueError('hsp, hpi, elp, point_names, fif, bvct must all ' - 'be None if egi is not None.') - _check_fname(egi, overwrite='read', must_exist=True) - - root = ElementTree.parse(egi).getroot() - ns = root.tag[root.tag.index('{'):root.tag.index('}') + 1] - sensors = root.find('%ssensorLayout/%ssensors' % (ns, ns)) - fids = dict() - dig_ch_pos = dict() - - fid_name_map = {'Nasion': 'nasion', - 'Right periauricular point': 'rpa', - 'Left periauricular point': 'lpa'} - - scale = dict(mm=1e-3, cm=1e-2, auto=1e-2, m=1) - if unit not in scale: - raise ValueError("Unit needs to be one of %s, not %r" % - (sorted(scale.keys()), unit)) - - for s in sensors: - name, number, kind = s[0].text, int(s[1].text), int(s[2].text) - coordinates = np.array([float(s[3].text), float(s[4].text), - float(s[5].text)]) - - coordinates *= scale[unit] - - # EEG Channels - if kind == 0: - dig_ch_pos['EEG %03d' % number] = coordinates - # Reference - elif kind == 1: - dig_ch_pos['EEG %03d' % - (len(dig_ch_pos.keys()) + 1)] = coordinates - # Fiducials - elif kind == 2: - fid_name = fid_name_map[name] - fids[fid_name] = coordinates - # Unknown - else: - warn('Unknown sensor type %s detected. Skipping sensor...' - 'Proceed with caution!' % kind) - - fids = [fids[key] for key in ('nasion', 'lpa', 'rpa')] - coord_frame = 'unknown' + data = _read_dig_montage_egi( + fname=egi, + _scaling=_get_scaling(unit, EGI_SCALE), + _all_data_kwargs_are_none=all( + x is None for x in (hsp, hpi, elp, point_names, fif, bvct)) + ) elif bvct is not None: - if not all(x is None for x in (hsp, hpi, elp, point_names, fif, egi)): - raise ValueError('hsp, hpi, elp, point_names, fif, egi must all ' - 'be None if bvct is not None.') - _check_fname(bvct, overwrite='read', must_exist=True) - - root = ElementTree.parse(bvct).getroot() - sensors = root.find('CapTrakElectrodeList') - - fids = {} - dig_ch_pos = {} - - fid_name_map = {'Nasion': 'nasion', 'RPA': 'rpa', 'LPA': 'lpa'} - - # CapTrak is natively in mm - scale = dict(mm=1e-3, cm=1e-2, auto=1e-3, m=1) - if unit not in scale: - raise ValueError("Unit needs to be one of %s, not %r" % - (sorted(scale.keys()), unit)) - if unit not in ['mm', 'auto']: - warn('Using "{}" as unit for BVCT file. BVCT files are usually ' - 'specified in "mm". This might lead to errors.'.format(unit), - RuntimeWarning) - - for s in sensors: - name = s.find('Name').text - - # Need to prune "GND" and "REF": these are not included in the raw - # data and will raise errors when we try to do raw.set_montage(...) - # XXX eventually this should be stored in ch['loc'][3:6] - # but we don't currently have such capabilities here - if name in ['GND', 'REF']: - continue - - fid = name in fid_name_map - coordinates = np.array([float(s.find('X').text), - float(s.find('Y').text), - float(s.find('Z').text)]) + data = _read_dig_montage_bvct( + fname=bvct, + unit=unit, # XXX: this should change + _all_data_kwargs_are_none=all( + x is None for x in (hsp, hpi, elp, point_names, fif, egi)) + ) - coordinates *= scale[unit] - - # Fiducials - if fid: - fid_name = fid_name_map[name] - fids[fid_name] = coordinates - # EEG Channels - else: - dig_ch_pos[name] = coordinates - - fids = [fids[key] for key in ('nasion', 'lpa', 'rpa')] - coord_frame = 'unknown' else: - fids = [None] * 3 - dig_ch_pos = None - scale = dict(mm=1e-3, cm=1e-2, auto=1e-3, m=1) - if unit not in scale: - raise ValueError("Unit needs to be one of %s, not %r" % - (sorted(scale.keys()), unit)) - + # XXX: This should also become a function + _scaling = _get_scaling(unit, NUMPY_DATA_SCALE), # HSP if isinstance(hsp, str): hsp = _read_dig_points(hsp, unit=unit) elif hsp is not None: - hsp *= scale[unit] + hsp *= _scaling # HPI if isinstance(hpi, str): @@ -826,18 +913,29 @@ def read_dig_montage(hsp=None, hpi=None, elp=None, point_names=None, # ELP if isinstance(elp, str): elp = _read_dig_points(elp, unit=unit) - elif elp is not None and scale[unit]: - elp *= scale[unit] - coord_frame = 'unknown' + elif elp is not None: + elp *= _scaling + + data = Bunch( + nasion=None, lpa=None, rpa=None, + hsp=hsp, elp=elp, coord_frame='unknown', + dig_ch_pos=None, hpi=hpi, point_names=point_names, + ) + + if any(x is None for x in (data.nasion, data.rpa, data.lpa)) and transform: + data = _fix_data_fiducials(data) - # Transform digitizer coordinates to neuromag space - out = DigMontage(hsp, hpi, elp, point_names, fids[0], fids[1], fids[2], - dig_ch_pos=dig_ch_pos, coord_frame=coord_frame) - if fif is None and transform: # only need to do this for non-Neuromag - out.transform_to_head() - if dev_head_t: - out.compute_dev_head_t() - return out + point_names = data.pop('point_names') + data['hpi_dev'] = data['hpi'] + data['hpi'] = data.pop('elp') + data['ch_pos'] = data.pop('dig_ch_pos') + montage = make_dig_montage( + **data, transform_to_head=transform, + compute_dev_head_t=dev_head_t, + ) + + montage._point_names = point_names # XXX: hack this should go!! + return montage def _set_montage(info, montage, update_ch_names=False, set_dig=True): @@ -929,20 +1027,25 @@ def _set_montage(info, montage, update_ch_names=False, set_dig=True): 'left untouched.') elif isinstance(montage, DigMontage): + if set_dig: - info['dig'] = montage._get_dig() + info['dig'] = montage.dig if montage.dev_head_t is not None: info['dev_head_t']['trans'] = montage.dev_head_t - if montage.dig_ch_pos is not None: # update channel positions, too - eeg_ref_pos = montage.dig_ch_pos.get('EEG000', np.zeros(3)) + if montage.ch_names: # update channel positions, too + dig_ch_pos = dict(zip(montage.ch_names, [ + d['r'] for d in montage.dig + if d['kind'] == FIFF.FIFFV_POINT_EEG + ])) + eeg_ref_pos = dig_ch_pos.get('EEG000', np.zeros(3)) did_set = np.zeros(len(info['ch_names']), bool) is_eeg = np.zeros(len(info['ch_names']), bool) is_eeg[pick_types(info, meg=False, eeg=True, exclude=())] = True - for ch_name, ch_pos in montage.dig_ch_pos.items(): - if ch_name == 'EEG000': + for ch_name, ch_pos in dig_ch_pos.items(): + if ch_name == 'EEG000': # what if eeg ref. has different name? continue if ch_name not in info['ch_names']: raise RuntimeError('Montage channel %s not found in info' @@ -958,6 +1061,7 @@ def _set_montage(info, montage, update_ch_names=False, set_dig=True): if len(did_not_set) > 0: warn('Did not set %s channel positions:\n%s' % (len(did_not_set), ', '.join(did_not_set))) + elif montage is None: for ch in info['chs']: ch['loc'] = np.full(12, np.nan) diff --git a/mne/channels/tests/test_montage.py b/mne/channels/tests/test_montage.py index 41a56a96c3f..745a32e1795 100644 --- a/mne/channels/tests/test_montage.py +++ b/mne/channels/tests/test_montage.py @@ -19,14 +19,16 @@ from mne import create_info, EvokedArray, read_evokeds, __file__ as _mne_file from mne.channels import (Montage, read_montage, read_dig_montage, - get_builtin_montages) + get_builtin_montages, DigMontage) from mne.channels.montage import _set_montage +from mne.channels._dig_montage_utils import _transform_to_head_call +from mne.channels._dig_montage_utils import _fix_data_fiducials from mne.utils import (_TempDir, run_tests_if_main, assert_dig_allclose, - object_diff) + object_diff, Bunch) from mne.bem import _fit_sphere -from mne.coreg import fit_matched_points from mne.transforms import apply_trans, get_ras_to_neuromag_trans from mne.io.constants import FIFF +from mne.digitization import Digitization from mne.digitization._utils import _read_dig_points from mne.viz._3d import _fiducial_coords @@ -376,40 +378,49 @@ def test_read_dig_montage(): elp_points = _read_dig_points(elp) hsp_points = _read_dig_points(hsp) hpi_points = read_mrk(hpi) - assert_equal(montage.point_names, names) - assert_array_equal(montage.elp, elp_points) - assert_array_equal(montage.hsp, hsp_points) - assert_array_equal(montage.hpi, hpi_points) + with pytest.deprecated_call(): + assert_equal(montage.point_names, names) + assert_array_equal(montage.elp, elp_points) + assert_array_equal(montage.hsp, hsp_points) assert (montage.dev_head_t is None) montage = read_dig_montage(hsp, hpi, elp, names, transform=True, dev_head_t=True) # check coordinate transformation # nasion - assert_almost_equal(montage.nasion[0], 0) - assert_almost_equal(montage.nasion[2], 0) + with pytest.deprecated_call(): + assert_almost_equal(montage.nasion[0], 0) + assert_almost_equal(montage.nasion[2], 0) # lpa and rpa - assert_allclose(montage.lpa[1:], 0, atol=1e-16) - assert_allclose(montage.rpa[1:], 0, atol=1e-16) + with pytest.deprecated_call(): + assert_allclose(montage.lpa[1:], 0, atol=1e-16) + assert_allclose(montage.rpa[1:], 0, atol=1e-16) # device head transform - dev_head_t = fit_matched_points(tgt_pts=montage.elp, - src_pts=montage.hpi, out='trans') - assert_array_equal(montage.dev_head_t, dev_head_t) + + EXPECTED_DEV_HEAD_T = np.array( + [[-3.72201691e-02, -9.98212167e-01, -4.67667497e-02, -7.31583414e-04], + [8.98064989e-01, -5.39382685e-02, 4.36543170e-01, 1.60134431e-02], + [-4.38285221e-01, -2.57513699e-02, 8.98466990e-01, 6.13035748e-02], + [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]] + ) + assert_allclose(montage.dev_head_t, EXPECTED_DEV_HEAD_T, atol=1e-7) # Digitizer as array m2 = read_dig_montage(hsp_points, hpi_points, elp_points, names, unit='m') - assert_array_equal(m2.hsp, montage.hsp) + with pytest.deprecated_call(): + assert_array_equal(m2.hsp, montage.hsp) m3 = read_dig_montage(hsp_points * 1000, hpi_points, elp_points * 1000, names) - assert_allclose(m3.hsp, montage.hsp) + with pytest.deprecated_call(): + assert_allclose(m3.hsp, montage.hsp) # test unit parameter and .mat support tempdir = _TempDir() mat_hsp = op.join(tempdir, 'test.mat') savemat(mat_hsp, dict(Points=(1000 * hsp_points).T), oned_as='row') montage_cm = read_dig_montage(mat_hsp, hpi, elp, names, unit='cm') - assert_allclose(montage_cm.hsp, montage.hsp * 10.) - assert_allclose(montage_cm.elp, montage.elp * 10.) - assert_array_equal(montage_cm.hpi, montage.hpi) + with pytest.deprecated_call(): + assert_allclose(montage_cm.hsp, montage.hsp * 10.) + assert_allclose(montage_cm.elp, montage.elp * 10.) pytest.raises(ValueError, read_dig_montage, hsp, hpi, elp, names, unit='km') # extra columns @@ -424,8 +435,9 @@ def test_read_dig_montage(): fout.write(line.rstrip() + b' 0.0 0.0 0.0\n') with pytest.warns(RuntimeWarning, match='Found .* columns instead of 3'): montage_extra = read_dig_montage(extra_hsp, hpi, elp, names) - assert_allclose(montage_extra.hsp, montage.hsp) - assert_allclose(montage_extra.elp, montage.elp) + with pytest.deprecated_call(): + assert_allclose(montage_extra.hsp, montage.hsp) + assert_allclose(montage_extra.elp, montage.elp) def test_set_dig_montage(): @@ -495,9 +507,6 @@ def test_fif_dig_montage(): raw_bv.add_channels([raw_bv_2]) for ii in range(2): - if ii == 1: - dig_montage.transform_to_head() # should have no meaningful effect - # Set the montage raw_bv.set_montage(dig_montage) @@ -534,14 +543,15 @@ def test_egi_dig_montage(): fname_temp = op.join(temp_dir, 'egi_test.fif') _check_roundtrip(dig_montage, fname_temp) - # Test coordinate transform - dig_montage.transform_to_head() - # nasion - assert_almost_equal(dig_montage.nasion[0], 0) - assert_almost_equal(dig_montage.nasion[2], 0) - # lpa and rpa - assert_allclose(dig_montage.lpa[1:], 0, atol=1e-16) - assert_allclose(dig_montage.rpa[1:], 0, atol=1e-16) + with pytest.deprecated_call(): + # Test coordinate transform + # dig_montage.transform_to_head() # XXX: this call had no effect!! + # nasion + assert_almost_equal(dig_montage.nasion[0], 0) + assert_almost_equal(dig_montage.nasion[2], 0) + # lpa and rpa + assert_allclose(dig_montage.lpa[1:], 0, atol=1e-16) + assert_allclose(dig_montage.rpa[1:], 0, atol=1e-16) # Test accuracy and embedding within raw object raw_egi = read_raw_egi(egi_raw_fname, channel_naming='EEG %03d') @@ -570,18 +580,17 @@ def test_bvct_dig_montage(): fname_temp = op.join(temp_dir, 'bvct_test.fif') _check_roundtrip(dig_montage, fname_temp) - # Test coordinate transform - dig_montage.transform_to_head() - # nasion - assert_almost_equal(dig_montage.nasion[0], 0) - assert_almost_equal(dig_montage.nasion[2], 0) - # lpa and rpa - assert_allclose(dig_montage.lpa[1:], 0, atol=1e-16) - assert_allclose(dig_montage.rpa[1:], 0, atol=1e-16) + with pytest.deprecated_call(): + # nasion + assert_almost_equal(dig_montage.nasion[0], 0) + assert_almost_equal(dig_montage.nasion[2], 0) + # lpa and rpa + assert_allclose(dig_montage.lpa[1:], 0, atol=1e-16) + assert_allclose(dig_montage.rpa[1:], 0, atol=1e-16) # Test accuracy and embedding within raw object raw_bv = read_raw_brainvision(bv_raw_fname) - with pytest.warns(RuntimeWarning, match='Did not set 3 channel pos'): + with pytest.warns(RuntimeWarning, match='Did not set.*channel pos'): raw_bv.set_montage(dig_montage) test_raw_bv = read_raw_fif(bv_fif_fname) @@ -613,15 +622,18 @@ def test_set_montage(): def _check_roundtrip(montage, fname): """Check roundtrip writing.""" - assert_equal(montage.coord_frame, 'head') + with pytest.deprecated_call(): + assert_equal(montage.coord_frame, 'head') montage.save(fname) montage_read = read_dig_montage(fif=fname) assert_equal(str(montage), str(montage_read)) - for kind in ('elp', 'hsp', 'nasion', 'lpa', 'rpa'): - if getattr(montage, kind) is not None: - assert_allclose(getattr(montage, kind), - getattr(montage_read, kind), err_msg=kind) - assert_equal(montage_read.coord_frame, 'head') + with pytest.deprecated_call(): + for kind in ('elp', 'hsp', 'nasion', 'lpa', 'rpa'): + if getattr(montage, kind, None) is not None: + assert_allclose(getattr(montage, kind), + getattr(montage_read, kind), err_msg=kind) + with pytest.deprecated_call(): + assert_equal(montage_read.coord_frame, 'head') def _fake_montage(ch_names): @@ -798,4 +810,57 @@ def test_setting_hydrocel_montage(): assert actual == expected +def test_dig_dev_head_t_regression(): + """Test deprecated compute_dev_head_t behavior.""" + def _read_dig_montage( + hsp=None, hpi=None, elp=None, point_names=None, unit='auto', + fif=None, egi=None, bvct=None, transform=True, dev_head_t=False, + ): + """Unfolds the `read_dig_montage` old behavior of the call below. + + montage = read_dig_montage(hsp, hpi, elp, names, + transform=True, dev_head_t=False) + """ + assert isinstance(hsp, str), 'original call hsp was string' + assert op.splitext(hpi)[-1] == '.sqd', 'original call hpi was .sqd' + assert isinstance(elp, str), 'original call elp was string' + + hsp = _read_dig_points(hsp, unit=unit) + hpi = read_mrk(hpi) + elp = _read_dig_points(elp, unit=unit) + + data = Bunch(nasion=None, lpa=None, rpa=None, + hsp=hsp, hpi=hpi, elp=elp, coord_frame='unknown', + point_names=point_names, dig_ch_pos=None) + + data = _fix_data_fiducials(data) + data = _transform_to_head_call(data) + with pytest.deprecated_call(): + montage = DigMontage(**data) + + return montage + + EXPECTED_DEV_HEAD_T = \ + [[-3.72201691e-02, -9.98212167e-01, -4.67667497e-02, -7.31583414e-04], + [8.98064989e-01, -5.39382685e-02, 4.36543170e-01, 1.60134431e-02], + [-4.38285221e-01, -2.57513699e-02, 8.98466990e-01, 6.13035748e-02], + [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]] + + names = ['nasion', 'lpa', 'rpa', '1', '2', '3', '4', '5'] + montage = _read_dig_montage( + hsp, hpi, elp, names, transform=True, dev_head_t=False) + + assert montage.dev_head_t is None + with pytest.deprecated_call(): + montage.compute_dev_head_t() + assert_allclose(montage.dev_head_t, EXPECTED_DEV_HEAD_T, atol=1e-7) + + +def test_make_dig_montage_errors(): + """Test proper error messaging.""" + with pytest.raises(ValueError, match='does not match the number'): + _ = DigMontage(ch_names=['foo', 'bar'], dig=Digitization()) + with pytest.raises(TypeError, match='must be an instance of Digitization'): + _ = DigMontage(ch_names=['foo', 'bar'], dig=None) + run_tests_if_main() diff --git a/mne/coreg.py b/mne/coreg.py index 97d9a89b846..daf3744f26d 100644 --- a/mne/coreg.py +++ b/mne/coreg.py @@ -289,6 +289,7 @@ def _trans_from_params(param_info, params): return trans +# XXX this function should be moved out of coreg as used elsewhere def fit_matched_points(src_pts, tgt_pts, rotate=True, translate=True, scale=False, tol=None, x0=None, out='trans', weights=None): diff --git a/mne/digitization/_utils.py b/mne/digitization/_utils.py index 316eded3486..d849c3c66e3 100644 --- a/mne/digitization/_utils.py +++ b/mne/digitization/_utils.py @@ -30,6 +30,7 @@ from ..transforms import combine_transforms from ..transforms import invert_transform from ..transforms import _to_const +from ..transforms import _str_to_frame from ..utils.check import _check_option from .. import __version__ @@ -182,9 +183,9 @@ def _write_dig_points(fname, dig_points): raise ValueError(msg) -# XXX: all points are supposed to be in FIFFV_COORD_HEAD def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, - extra_points=None, dig_ch_pos=None): + extra_points=None, dig_ch_pos=None, + coord_frame='head'): """Construct digitizer info for the info. Parameters @@ -201,12 +202,21 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, Points designed as the headshape points. dig_ch_pos : dict Dict of EEG channel positions. + coord_frame : str + The coordinate frame of the points. Usually this is "unknown" + for native digitizer space. Defaults to "head". Returns ------- dig : list of dicts A container of DigPoints to be added to the info['dig']. """ + if not isinstance(coord_frame, str) or coord_frame not in _str_to_frame: + raise ValueError('coord_frame must be one of %s, got %s' + % (sorted(_str_to_frame.keys()), coord_frame)) + else: + coord_frame = _str_to_frame[coord_frame] + dig = [] if lpa is not None: lpa = np.asarray(lpa) @@ -215,7 +225,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, % (lpa.shape,)) dig.append({'r': lpa, 'ident': FIFF.FIFFV_POINT_LPA, 'kind': FIFF.FIFFV_POINT_CARDINAL, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) if nasion is not None: nasion = np.asarray(nasion) if nasion.shape != (3,): @@ -223,7 +233,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, % (nasion.shape,)) dig.append({'r': nasion, 'ident': FIFF.FIFFV_POINT_NASION, 'kind': FIFF.FIFFV_POINT_CARDINAL, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) if rpa is not None: rpa = np.asarray(rpa) if rpa.shape != (3,): @@ -231,7 +241,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, % (rpa.shape,)) dig.append({'r': rpa, 'ident': FIFF.FIFFV_POINT_RPA, 'kind': FIFF.FIFFV_POINT_CARDINAL, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) if hpi is not None: hpi = np.asarray(hpi) if hpi.ndim != 2 or hpi.shape[1] != 3: @@ -240,7 +250,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, for idx, point in enumerate(hpi): dig.append({'r': point, 'ident': idx + 1, 'kind': FIFF.FIFFV_POINT_HPI, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) if extra_points is not None: extra_points = np.asarray(extra_points) if extra_points.shape[1] != 3: @@ -249,7 +259,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, for idx, point in enumerate(extra_points): dig.append({'r': point, 'ident': idx + 1, 'kind': FIFF.FIFFV_POINT_EXTRA, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) if dig_ch_pos is not None: keys = sorted(dig_ch_pos.keys()) try: # use the last 3 as int if possible (e.g., EEG001->1) @@ -263,7 +273,7 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None, for key, ident in zip(keys, idents): dig.append({'r': dig_ch_pos[key], 'ident': ident, 'kind': FIFF.FIFFV_POINT_EEG, - 'coord_frame': FIFF.FIFFV_COORD_HEAD}) + 'coord_frame': coord_frame}) return _format_dig_points(dig) diff --git a/mne/gui/_file_traits.py b/mne/gui/_file_traits.py index f99ba3c11a8..e12ce444b5e 100644 --- a/mne/gui/_file_traits.py +++ b/mne/gui/_file_traits.py @@ -300,20 +300,12 @@ def _get__info(self): if len(dir_tree_find(tree, FIFF.FIFFB_MEAS_INFO)) > 0: info = read_info(self.file, verbose=False) elif len(dir_tree_find(tree, FIFF.FIFFB_ISOTRAK)) > 0: - info = read_dig_montage(fif=self.file) + info = read_dig_montage(fif=self.file, transform=True) if isinstance(info, DigMontage): - info.transform_to_head() - digs = list() - _append_fiducials(digs, info.lpa, info.nasion, info.rpa) - for idx, pos in enumerate(info.hsp): - dig = {'coord_frame': FIFF.FIFFV_COORD_HEAD, - 'ident': idx, - 'kind': FIFF.FIFFV_POINT_EXTRA, - 'r': pos} - digs.append(dig) + dig = info.dig info = _empty_info(1) - info['dig'] = digs + info['dig'] = dig elif info is None or info['dig'] is None: error(None, "The selected FIFF file does not contain " "digitizer information. Please select a different " diff --git a/mne/io/array/tests/test_array.py b/mne/io/array/tests/test_array.py index b04b68c98c8..d249b1b1e67 100644 --- a/mne/io/array/tests/test_array.py +++ b/mne/io/array/tests/test_array.py @@ -10,12 +10,13 @@ import pytest import matplotlib.pyplot as plt -from mne import find_events, Epochs, pick_types, channels +from mne import find_events, Epochs, pick_types from mne.io import read_raw_fif from mne.io.array import RawArray from mne.io.tests.test_raw import _test_raw_reader from mne.io.meas_info import create_info, _kind_dict from mne.utils import requires_version, run_tests_if_main +from mne.channels import make_dig_montage base_dir = op.join(op.dirname(__file__), '..', '..', 'tests', 'data') fif_fname = op.join(base_dir, 'test_raw.fif') @@ -158,15 +159,14 @@ def test_array_raw(): n_elec = 10 ts_size = 10000 Fs = 512. - elec_labels = [str(i) for i in range(n_elec)] - elec_coords = np.random.randint(60, size=(n_elec, 3)).tolist() + ch_names = [str(i) for i in range(n_elec)] + ch_pos_loc = np.random.randint(60, size=(n_elec, 3)).tolist() - electrode = np.random.rand(n_elec, ts_size) - dig_ch_pos = dict(zip(elec_labels, elec_coords)) - mon = channels.DigMontage(dig_ch_pos=dig_ch_pos) - info = create_info(elec_labels, Fs, 'ecog', montage=mon) + data = np.random.rand(n_elec, ts_size) + montage = make_dig_montage(ch_pos=dict(zip(ch_names, ch_pos_loc))) + info = create_info(ch_names, Fs, 'ecog', montage=montage) - raw = RawArray(electrode, info) + raw = RawArray(data, info) raw.plot_psd(average=False) # looking for inexistent layout raw.plot_psd_topo() diff --git a/mne/io/fieldtrip/utils.py b/mne/io/fieldtrip/utils.py index dc5863c0344..371b33bb9c7 100644 --- a/mne/io/fieldtrip/utils.py +++ b/mne/io/fieldtrip/utils.py @@ -7,7 +7,7 @@ from ..meas_info import create_info from ...transforms import rotation3d_align_z_axis -from ...channels import DigMontage +from ...channels import make_dig_montage from ..constants import FIFF from ...utils import warn, _check_pandas_installed from ..pick import pick_info @@ -148,8 +148,9 @@ def _create_montage(ft_struct): if (len(montage_ch_names) > 0 and len(montage_pos) > 0 and len(montage_ch_names) == len(montage_pos)): - montage = DigMontage( - dig_ch_pos=dict(zip(montage_ch_names, montage_pos))) + montage = make_dig_montage( + ch_pos=dict(zip(montage_ch_names, montage_pos)) + ) return montage diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py index 0db3a1c6646..4417f255da0 100644 --- a/mne/viz/_3d.py +++ b/mne/viz/_3d.py @@ -2584,15 +2584,15 @@ def snapshot_brain_montage(fig, montage, hide_sensors=True): im : array, shape (m, n, 3) The screenshot of the current scene view """ - from ..channels import Montage, DigMontage + from ..channels import DigMontage from .. import Info # Update the backend from .backends.renderer import _Renderer if fig is None: raise ValueError('The figure must have a scene') - if isinstance(montage, (Montage, DigMontage)): - chs = montage.dig_ch_pos + if isinstance(montage, DigMontage): + chs = montage._get_ch_pos() ch_names, xyz = zip(*[(ich, ixyz) for ich, ixyz in chs.items()]) elif isinstance(montage, Info): xyz = [ich['loc'][:3] for ich in montage['chs']] diff --git a/mne/viz/montage.py b/mne/viz/montage.py index 13f55ee573e..0a1617e01d9 100644 --- a/mne/viz/montage.py +++ b/mne/viz/montage.py @@ -35,7 +35,7 @@ def plot_montage(montage, scale_factor=20, show_names=True, kind='topomap', ch_names = montage.ch_names title = montage.kind elif isinstance(montage, DigMontage): - ch_names = montage.point_names + ch_names = montage._point_names title = None else: raise TypeError("montage must be an instance of " diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 01bbb30af9e..4ea57752baa 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -28,15 +28,14 @@ mat = loadmat(mne.datasets.misc.data_path() + '/ecog/sample_ecog.mat') ch_names = mat['ch_names'].tolist() elec = mat['elec'] # electrode positions given in meters -dig_ch_pos = dict(zip(ch_names, elec)) -mon = mne.channels.DigMontage(dig_ch_pos=dig_ch_pos) +montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec))) print('Created %s channel positions' % len(ch_names)) ############################################################################### # Now that we have our electrode positions in MRI coordinates, we can create # our measurement info structure. -info = mne.create_info(ch_names, 1000., 'ecog', montage=mon) +info = mne.create_info(ch_names, 1000., 'ecog', montage=montage) ############################################################################### # We can then plot the locations of our electrodes on our subject's brain. @@ -59,7 +58,7 @@ fig_scatter = plot_alignment(info, subject='sample', subjects_dir=subjects_dir, surfaces='pial') mne.viz.set_3d_view(fig_scatter, 200, 70) -xy, im = snapshot_brain_montage(fig_scatter, mon) +xy, im = snapshot_brain_montage(fig_scatter, montage) # Convert from a dictionary to array to plot xy_pts = np.vstack([xy[ch] for ch in info['ch_names']])