MNT: Change DigMontage representation to use Digitization (mne-tools#…

…6639)

* WIP: make transformations private [run CIs]

* cosmit and notes

* WIP: move transforms 

This commit corresponds to this PR:
massich#27

* WIP: move the transformations logic into free helper functions (mne-tools#28)

* WIP: Add dig representation

* WIP: Remove transforms from  __init__

See massich#29 for details

* FIX dictionary unpacking

* fix rebase

* WIP: Keep dig and chnames only 

from massich#30

* fix 052193c

* use deprecated property

* wip: make point_names private + property

* wip: remove deprecation warnings

* wip: add ch_names attribute

* wip: do not use point_names in __repr__

* wip: remove _get_dig()

* wip with alex!!!!

* wip: fix wip

* fix pep

* xxxx:

* wip: deprecate compute_dev_head_t

* XXX: update XXX comments

* wip: add meaningful deprecation message

* wip: use new DigMontage constructor

* WIP-TST: deprecated DigMontage contruction

* TST: use atol for dev_head_t testing

* WIP: ADD make_dig_montage (and make it great again)

* wip

* fix tests

* fix test?

* don't copy

* cleanup

* more fixes

* fix

* fix fieldtrip

* FIX: Link

* update whatsnew and trigger CIs

* use master whatsnew

* update the new whatsnew

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`_

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',

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

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']])

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)

mne/channels/_dig_montage_utils.py

@@ -0,0 +1,306 @@
+# Authors: Alexandre Gramfort <[email protected]>
+#          Denis Engemann <[email protected]>
+#          Martin Luessi <[email protected]>
+#          Eric Larson <[email protected]>
+#          Marijn van Vliet <[email protected]>
+#          Jona Sassenhagen <[email protected]>
+#          Teon Brooks <[email protected]>
+#          Christian Brodbeck <[email protected]>
+#          Stefan Appelhoff <[email protected]>
+#          Joan Massich <[email protected]>
+#
+# 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,
+    )