warp.py

from copy import deepcopy

import mne
from mne.bem import _fit_sph_harm_to_headshape, _transform_sph
from mne.surface import _tessellate_sphere_surf

raw = mne.io.read_raw_fif(mne.datasets.sample.data_path() +
                          '/MEG/sample/sample_audvis_raw.fif')

# For a visual explanation of the order, see:
#     https://en.wikipedia.org/wiki/Spherical_harmonics
# Rows correspond to orders, where 0 is just a sphere, i.e. just a
# radius fit.

order = 4  # 0 = translation + radius
_, origin, _, sph_coeffs, dig = _fit_sph_harm_to_headshape(
    raw.info, ['extra', 'eeg'], order=order, verbose=True)

print('Loading head (dense) and sphere surfaces in head coordinates')
data_path = mne.datasets.sample.data_path()
head = mne.read_bem_surfaces(
    data_path + '/subjects/sample/bem/sample-head.fif')[0]
trans = mne.read_trans(data_path + '/MEG/sample/sample_audvis_raw-trans.fif')
mne.transform_surface_to(head, 'head', trans)
sphere = _tessellate_sphere_surf(4)

print('Transforming surfaces using spherical harmonics')
head_trans = deepcopy(head)
head_trans['rr'] = _transform_sph(head['rr'], sph_coeffs, origin)
sphere['rr'] = _transform_sph(sphere['rr'], sph_coeffs, origin)

print('Plot dig (black), original (gray), transform (purple), & fit (green)')
from mayavi import mlab  # noqa
fig = mlab.figure(size=(800, 600), bgcolor=(1., 1., 1.))
for surf, color, opacity in zip((head_trans, head, sphere),
                                ((0.3, 1., 0.5), (0.3, 0.3, 0.3), (1, 0., 1)),
                                (0.8, 0.5, 0.3)):
    mesh = mlab.pipeline.triangular_mesh_source(
        *surf['rr'].T, triangles=surf['tris'])
    mesh.data.point_data.normals = surf['nn']
    mesh.data.cell_data.normals = None
    surf = mlab.pipeline.surface(mesh, figure=fig, reset_zoom=True,
                                 opacity=opacity, color=color)
    surf.actor.property.backface_culling = True
mlab.points3d(dig[:, 0], dig[:, 1], dig[:, 2], color=(0.1, 0.1, 0.1),
              scale_factor=0.005)
mlab.view(45, 90)