head_trans.py

# -*- coding: utf-8 -*-

# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD (3-clause)

from __future__ import print_function

from copy import deepcopy
import time

import numpy as np

import mne
from mne.preprocessing import maxwell_filter
from mne.simulation import simulate_raw
from mne.transforms import transform_surface_to, Transform, read_trans

data_path = mne.datasets.sample.data_path()

simulate = False # re-simulate data?

n_sim = 8  # 1 stationary plus 7 rots
origin = (0., 0., 0.05)
translation = (0., 0., 0.04)  # down the helmet just a bit
raw_fnames = ['sim_%02d_raw.fif' % ri for ri in range(n_sim)]

trans_fname = data_path + '/MEG/sample/sample_audvis_raw-trans.fif'
subjects_dir = data_path + '/subjects'
bem_fname = subjects_dir + '/sample/bem/sample-5120-bem-sol.fif'
src_fname = subjects_dir + '/sample/bem/sample-oct-6-src.fif'

src = mne.read_source_spaces(src_fname)
label = mne.read_labels_from_annot('sample', 'aparc.a2009s', 'lh',
                                   regexp='G_temp_sup-Plan_tempo',
                                   subjects_dir=subjects_dir)
assert len(label) == 1
label = label[0]
vertices = [np.intersect1d(label.vertices, src[0]['vertno'])[2:3],
            np.array([], int)]
trans = read_trans(trans_fname)
src_head = transform_surface_to(deepcopy(src[0]), 'head', trans)
pos = src_head['rr'][vertices[0][0]]
ori = src_head['nn'][vertices[0][0]]
amp = 10e-9  # nAm
if simulate:
    raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'

    # let's make rotation matrices about each axis, plus some compound ones
    phi = np.deg2rad(30)
    x_rot = np.array([[1, 0, 0],
                      [0, np.cos(phi), -np.sin(phi)],
                      [0, np.sin(phi), np.cos(phi)]])
    y_rot = np.array([[np.cos(phi), 0, np.sin(phi)],
                      [0, 1, 0],
                      [-np.sin(phi), 0, np.cos(phi)]])
    z_rot = np.array([[np.cos(phi), -np.sin(phi), 0],
                      [np.sin(phi), np.cos(phi), 0],
                      [0, 0, 1]])
    xz_rot = np.dot(x_rot, z_rot)
    xmz_rot = np.dot(x_rot, z_rot.T)
    yz_rot = np.dot(y_rot, z_rot)
    mymz_rot = np.dot(y_rot.T, z_rot.T)

    # Create different head rotations, one per second
    rots = [x_rot, y_rot, z_rot, xz_rot, xmz_rot, yz_rot, mymz_rot]
    # The transpose of a rotation matrix is a rotation in the opposite dir
    rots = [np.eye(3)] + rots  # first one is in the destination position!
    assert len(rots) == n_sim

    raw = mne.io.Raw(raw_fname).crop(0, 10)
    raw.load_data().pick_types(meg=True, stim=True, copy=False)
    raw.add_proj([], remove_existing=True)

    # Let's activate a left auditory source
    data = np.zeros([sum(len(v) for v in vertices), int(raw.info['sfreq'])])
    activation = np.hanning(int(raw.info['sfreq'] * 0.2)) * amp
    t_offset = int(np.ceil(0.2 * raw.info['sfreq']))  # 200 ms in (after bl)
    data[:, t_offset:t_offset + len(activation)] = activation
    stc = mne.SourceEstimate(data, vertices, tmin=-0.2,
                             tstep=1. / raw.info['sfreq'])

    # Simulate the movement
    raw.info['dev_head_t']['trans'][:3, 3] = translation
    for ri, rot in enumerate(rots):
        raw.info['dev_head_t']['trans'][:3, :3] = rot
        raw_sim = simulate_raw(raw, stc, trans, src, bem_fname,
                               n_jobs=-1, verbose=True)
        raw_sim.save(raw_fnames[ri], buffer_size_sec=None, overwrite=True)
    stc.save('simulated_activation')
else:
    stc = mne.read_source_estimate('simulated_activation')
raw_fnames, stat_fname = raw_fnames[:-1], raw_fnames[-1]


def map_mn(raw, destination, origin):
    """Map raw data to a new dev_head_t"""
    from mne.forward._field_interpolation import _map_meg_channels
    from mne.io.pick import pick_info, pick_types
    picks_from = pick_types(raw.info, exclude='bads')
    picks_to = pick_types(raw.info, exclude=())
    info_from_meg = pick_info(raw.info, picks_from)
    info_to_meg = pick_info(raw.info, picks_to)
    info_to_meg['dev_head_t'] = destination
    m = _map_meg_channels(info_from_meg, info_to_meg, 'accurate', origin)
    raw = raw.copy().load_data()
    raw.info['dev_head_t'] = destination
    raw._data[picks_to] = np.dot(m, raw._data[picks_from])
    return raw

raws = [mne.io.Raw(r) for r in raw_fnames[1:]]  # exclude dest pos one

# stationary (at destination) version
raws_stat = [mne.io.Raw(raw_fnames[0])]

# transform others to the common destination (should have np.eye(3) as rot)
dest = Transform('meg', 'head', np.eye(4))
dest['trans'][:3, 3] = translation

print('Transforming with maxwell_filter ... ', end='')
t0 = time.time()
raws_mf = [maxwell_filter(r, origin=origin, destination=dest) for r in raws]
print(round(time.time() - t0))

print('Transforming with minimum norm ... ', end='')
t0 = time.time()
raws_mn = [map_mn(r, destination=dest, origin=origin) for r in raws]
print(round(time.time() - t0))

for use_raws, title in ((raws_stat, 'stationary'),
                        (raws, 'unprocessed'),
                        (raws_mf, 'maxwell_filter'),
                        (raws_mn, 'minimum norm'),
                        ):
    raw = mne.concatenate_raws(use_raws)
    events = mne.find_events(raw, 'STI 014')
    if title == 'stationary':
        assert len(events) == 10
    else:
        assert len(events) == 10 * (len(raw_fnames) - 1)
    epochs = mne.Epochs(raw, events)
    # "shrunk" breaks dipole fitting for some conditions
    method = 'shrunk' if title == 'minimum norm' else 'empirical'
    cov = mne.compute_covariance(epochs, tmax=0., method=method)
    evoked = epochs.average()
    # visualize
    evoked.plot_topomap(times=[0.0, 0.1, 0.2], title=title,
                        vmin=-200, vmax=200)
    # quantify
    dip = mne.fit_dipole(evoked.copy().crop(0.1, 0.1),
                         cov, bem_fname, trans=trans_fname)[0]
    dist = 1000 * np.sqrt(np.sum((dip.pos[0] - pos) ** 2))
    angle = np.rad2deg(np.arccos(np.dot(ori, dip.ori[0])))
    factor = 100 * (1. - dip.amplitude[0] / amp)
    print('Errors %s:' % title)
    print('  pos: %4.1f mm' % (dist,))
    print('  ori: %4.1f deg' % (angle,))
    print('  amp: %4.1f%%' % (factor,))