artemis_delaunay.py

"""
Create ARTEMIS123 helmet from sensor locations.
"""

import os.path as op
import numpy as np
from scipy.spatial import ConvexHull, Delaunay
import matplotlib.pyplot as plt

import mne
from mne.io.constants import FIFF

use_delaunay = True

art = mne.io.read_raw_artemis123(
    mne.datasets.testing.data_path() + '/ARTEMIS123/'
    'Artemis_Data_2017-04-04-15h-44m-22s_Motion_Translation-z.bin')
picks = mne.pick_types(art.info, meg=True, ref_meg=False)
picks = np.setdiff1d(picks, [119])
rr = np.array([c['loc'][:3] for ci, c in enumerate(art.info['chs'])
               if ci in picks])
if use_delaunay:
    rr_ch = rr.copy()
    # reduce to 2D
    rr_ch -= rr_ch.mean(axis=0, keepdims=True)
    assert (np.linalg.norm(rr_ch, axis=-1) < 0.12).all()
    assert (np.linalg.norm(rr_ch, axis=-1) > 0.05).all()
    rr_ch[:, :2] *= (rr_ch[:, 2].max() - rr_ch[:, [2]] + 0.05)
    plt.scatter(*rr_ch[:, :2].T)
    d = Delaunay(rr_ch[:, :2])
    assert d.simplices.shape[1] == 3
    del rr_ch
    tris = d.simplices
else:
    c = ConvexHull(rr)
    tris = c.simplices
    assert np.array_equal(np.unique(tris.ravel()), np.arange(118))
surf = dict(rr=rr, tris=tris, id=FIFF.FIFFV_MNE_SURF_MEG_HELMET,
            coord_frame=FIFF.FIFFV_COORD_DEVICE, sigma=1.0, np=len(rr),
            ntri=len(tris))
mne.write_bem_surfaces(op.join(op.dirname(mne.__file__), 'data', 'helmets',
                               'ARTEMIS123.fif.gz'), [surf])
mne.viz.plot_alignment(art.info, surfaces=(), verbose=True)