MRG: Add subject warping (mne-tools#3613)

* WIP: Add spherical harmonic transformations

MANIFEST.in

@@ -6,6 +6,7 @@ recursive-include examples *.txt
 recursive-include mne *.py
 recursive-include mne/data *.dat
 recursive-include mne/data *.sel
+recursive-include mne/data *.fif
 recursive-include mne/data *.fif.gz
 recursive-include mne/channels/data/montages *.elc
 recursive-include mne/channels/data/montages *.txt

doc/python_reference.rst

@@ -651,6 +651,8 @@ Functions:
    :toctree: generated/
    :template: function.rst
 
+   fit_sphere_to_headshape
+   get_fitting_dig
    make_watershed_bem
    make_flash_bem
    convert_flash_mris
@@ -664,6 +666,22 @@ Functions:
    restrict_forward_to_label
    restrict_forward_to_stc
 
+.. currentmodule:: mne.transforms
+
+.. autosummary::
+   :toctree: generated/
+   :template: class.rst
+
+   SphericalSurfaceWarp
+
+.. currentmodule:: mne.surface
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   complete_surface_info
+
 Inverse Solutions
 =================
 

doc/whats_new.rst

@@ -17,6 +17,8 @@ Changelog
 
     - Add filter plotting functions :func:`mne.viz.plot_filter` and :func:`mne.viz.plot_ideal_filter` as well as filter creation function :func:`mne.filter.create_filter` by `Eric Larson`_
 
+    - Add 3D thin-plate spline warping with spherical harmonic surface approximations in :class:`mne.transforms.SphericalSurfaceWarp` by `Eric Larson`_
+
 BUG
 ~~~
 

examples/visualization/plot_evoked_whitening.py

@@ -7,13 +7,13 @@
 matrix. It's an excellent quality check to see if baseline signals match
 the assumption of Gaussian white noise from which we expect values around
 0 with less than 2 standard deviations. Covariance estimation and diagnostic
-plots are based on [1].
+plots are based on [1]_.
 
 References
 ----------
-[1] Engemann D. and Gramfort A. (2015) Automated model selection in covariance
-    estimation and spatial whitening of MEG and EEG signals, vol. 108,
-    328-342, NeuroImage.
+.. [1] Engemann D. and Gramfort A. (2015) Automated model selection in
+    covariance estimation and spatial whitening of MEG and EEG signals, vol.
+    108, 328-342, NeuroImage.
 
 """
 # Authors: Alexandre Gramfort <[email protected]>

examples/visualization/plot_warp_surfaces.py

@@ -0,0 +1,83 @@
+"""
+=================================
+Warp head shapes between subjects
+=================================
+
+In this example, we warp data from one subject (sample) to another (fsaverage)
+using a spherical harmonic approximation of surfaces, followed by thin-plate
+spline (TPS) warping of the surface coordinates, as described in [1]_.
+
+References
+----------
+.. [1] Darvas F, Ermer JJ, Mosher JC, Leahy RM (2006). "Generic head models for
+       atlas-based EEG source analysis." Human Brain Mapping 27:129-143
+"""
+# Author: Eric Larson <[email protected]>
+#
+# License: BSD (3-clause)
+
+import os.path as op
+
+import mne
+
+fsaverage_path = op.join(op.dirname(mne.__file__), 'data', 'fsaverage')
+data_path = mne.datasets.sample.data_path()
+subjects_dir = op.join(data_path, 'subjects')
+
+###############################################################################
+# Load the source digitization, destination surfaces, and source surfaces
+# (with everything in head coordinates):
+
+# Digitization
+info = mne.io.read_info(op.join(data_path, 'MEG', 'sample',
+                                'sample_audvis_raw.fif'))
+hsp = mne.bem.get_fitting_dig(info, ('cardinal', 'extra'))
+
+# Destination head surface
+fsaverage_surfs = [mne.read_bem_surfaces(op.join(fsaverage_path,
+                                                 'fsaverage-%s.fif' % kind))[0]
+                   for kind in ('head', 'inner_skull-bem')]
+fsaverage_trans = mne.read_trans(op.join(fsaverage_path,
+                                         'fsaverage-trans.fif'))
+for surf in fsaverage_surfs:
+    mne.surface.transform_surface_to(surf, 'head', fsaverage_trans, copy=False)
+
+# Some source surfaces to transform as examples
+sample_surfs = mne.read_bem_surfaces(
+    op.join(subjects_dir, 'sample', 'bem', 'sample-5120-bem.fif'))
+sample_trans = mne.read_trans(op.join(data_path, 'MEG', 'sample',
+                                      'sample_audvis_raw-trans.fif'))
+for surf in sample_surfs:
+    mne.surface.transform_surface_to(surf, 'head', sample_trans, copy=False)
+
+###############################################################################
+# Transform surfaces using TPS warping:
+
+warp = mne.transforms.SphericalSurfaceWarp()
+warp.fit(source=hsp, destination=fsaverage_surfs[0]['rr'])
+for surf in sample_surfs:
+    surf['rr'] = warp.transform(surf['rr'])
+    # recompute our normals (only used for viz here)
+    mne.surface.complete_surface_info(surf, copy=False)
+hsp = warp.transform(hsp)
+
+###############################################################################
+# Plot the transformed surfaces and digitization (blue) on template (black):
+
+from mayavi import mlab  # noqa
+
+t_color = (0.1, 0.3, 1)
+fig = mlab.figure(size=(400, 600), bgcolor=(1., 1., 1.))
+for surf, color in zip(fsaverage_surfs + sample_surfs,
+                       [(0., 0., 0.)] * len(fsaverage_surfs) +
+                       [t_color] * len(sample_surfs)):
+    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=0.33, color=color)
+    surf.actor.property.backface_culling = True
+mlab.points3d(hsp[:, 0], hsp[:, 1], hsp[:, 2], color=t_color,
+              scale_factor=0.005, opacity=0.25)
+mlab.view(45, 90)

mne/bem.py

@@ -15,7 +15,6 @@
 import numpy as np
 from scipy import linalg
 
-from .utils import verbose, logger, run_subprocess, get_subjects_dir, warn
 from .transforms import _ensure_trans, apply_trans
 from .io import Info
 from .io.constants import FIFF
@@ -25,15 +24,13 @@
 from .io.tag import find_tag
 from .io.tree import dir_tree_find
 from .io.open import fiff_open
+from .utils import verbose, logger, run_subprocess, get_subjects_dir, warn
 from .externals.six import string_types
 
 
 # ############################################################################
 # Compute BEM solution
 
-# define VEC_DIFF(from,to,diff) {\
-# (diff)[X] = (to)[X] - (from)[X];\
-
 # The following approach is based on:
 #
 # de Munck JC: "A linear discretization of the volume conductor boundary
@@ -249,9 +246,9 @@ def _fwd_bem_ip_modify_solution(solution, ip_solution, ip_mult, n_tri):
 def _fwd_bem_linear_collocation_solution(m):
     """Compute the linear collocation potential solution"""
     # first, add surface geometries
-    from .surface import _complete_surface_info
+    from .surface import complete_surface_info
     for surf in m['surfs']:
-        _complete_surface_info(surf, verbose=False)
+        complete_surface_info(surf, copy=False, verbose=False)
 
     logger.info('Computing the linear collocation solution...')
     logger.info('    Matrix coefficients...')
@@ -801,7 +798,7 @@ def make_sphere_model(r0=(0., 0., 0.04), head_radius=0.09, info=None,
 
 
 # #############################################################################
-# Helpers
+# Sphere fitting
 
 _dig_kind_dict = {
     'cardinal': FIFF.FIFFV_POINT_CARDINAL,
@@ -825,7 +822,8 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
         Kind of digitization points to use in the fitting. These can be any
         combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
         be 'auto' (default), which will use only the 'extra' points if
-        enough are available, and if not, uses 'extra' and 'eeg' points.
+        enough (more than 10) are available, and if not, uses 'extra' and
+        'eeg' points.
     units : str
         Can be "m" (default) or "mm".
 
@@ -850,6 +848,44 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
     """
     if not isinstance(units, string_types) or units not in ('m', 'mm'):
         raise ValueError('units must be a "m" or "mm"')
+    radius, origin_head, origin_device = _fit_sphere_to_headshape(
+        info, dig_kinds)
+    if units == 'mm':
+        radius *= 1e3
+        origin_head *= 1e3
+        origin_device *= 1e3
+    return radius, origin_head, origin_device
+
+
+@verbose
+def get_fitting_dig(info, dig_kinds='auto', verbose=None):
+    """Get digitization points suitable for sphere fitting.
+
+    Parameters
+    ----------
+    info : instance of Info
+        The measurement info.
+    dig_kinds : list of str | str
+        Kind of digitization points to use in the fitting. These can be any
+        combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
+        be 'auto' (default), which will use only the 'extra' points if
+        enough (more than 10) are available, and if not, uses 'extra' and
+        'eeg' points.
+    verbose : bool, str or None
+        If not None, override default verbose level
+
+    Returns
+    -------
+    dig : array, shape (n_pts, 3)
+        The digitization points (in head coordinates) to use for fitting.
+
+    Notes
+    -----
+    This will exclude digitization locations that have ``z < 0 and y > 0``,
+    i.e. points on the nose and below the nose on the face.
+
+    .. versionadded:: 0.14
+    """
     if not isinstance(info, Info):
         raise TypeError('info must be an instance of Info not %s' % type(info))
     if info['dig'] is None:
@@ -859,10 +895,10 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
         if dig_kinds == 'auto':
             # try "extra" first
             try:
-                return fit_sphere_to_headshape(info, 'extra', units=units)
+                return get_fitting_dig(info, 'extra')
             except ValueError:
                 pass
-            return fit_sphere_to_headshape(info, ('extra', 'eeg'), units=units)
+            return get_fitting_dig(info, ('extra', 'eeg'))
         else:
             dig_kinds = (dig_kinds,)
     # convert string args to ints (first make dig_kinds mutable in case tuple)
@@ -880,7 +916,7 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
                            'contact mne-python developers')
 
     # exclude some frontal points (nose etc.)
-    hsp = [p for p in hsp if not (p[2] < 0 and p[1] > 0)]
+    hsp = np.array([p for p in hsp if not (p[2] < -1e-6 and p[1] > 1e-6)])
 
     if len(hsp) <= 10:
         kinds_str = ', '.join(['"%s"' % _dig_kind_rev[d]
@@ -891,40 +927,38 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
             raise ValueError(msg + ', at least 4 required')
         else:
             warn(msg + ', fitting may be inaccurate')
+    return hsp
+
 
+@verbose
+def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
+    """Fit a sphere to the given head shape."""
+    hsp = get_fitting_dig(info, dig_kinds)
     radius, origin_head = _fit_sphere(np.array(hsp), disp=False)
     # compute origin in device coordinates
     head_to_dev = _ensure_trans(info['dev_head_t'], 'head', 'meg')
     origin_device = apply_trans(head_to_dev, origin_head)
-    radius *= 1e3
-    origin_head *= 1e3
-    origin_device *= 1e3
-
-    logger.info('Fitted sphere radius:'.ljust(30) + '%0.1f mm' % radius)
+    logger.info('Fitted sphere radius:'.ljust(30) + '%0.1f mm'
+                % (radius * 1e3,))
     # 99th percentile on Wikipedia for Giabella to back of head is 21.7cm,
     # i.e. 108mm "radius", so let's go with 110mm
     # en.wikipedia.org/wiki/Human_head#/media/File:HeadAnthropometry.JPG
-    if radius > 110.:
+    if radius > 0.110:
         warn('Estimated head size (%0.1f mm) exceeded 99th '
-             'percentile for adult head size' % (radius,))
+             'percentile for adult head size' % (1e3 * radius,))
     # > 2 cm away from head center in X or Y is strange
-    if np.sqrt(np.sum(origin_head[:2] ** 2)) > 20:
+    if np.sqrt(np.sum(origin_head[:2] ** 2)) > 0.02:
         warn('(X, Y) fit (%0.1f, %0.1f) more than 20 mm from '
-             'head frame origin' % tuple(origin_head[:2]))
+             'head frame origin' % tuple(1e3 * origin_head[:2]))
     logger.info('Origin head coordinates:'.ljust(30) +
-                '%0.1f %0.1f %0.1f mm' % tuple(origin_head))
+                '%0.1f %0.1f %0.1f mm' % tuple(1e3 * origin_head))
     logger.info('Origin device coordinates:'.ljust(30) +
-                '%0.1f %0.1f %0.1f mm' % tuple(origin_device))
-    if units == 'm':
-        radius /= 1e3
-        origin_head /= 1e3
-        origin_device /= 1e3
-
+                '%0.1f %0.1f %0.1f mm' % tuple(1e3 * origin_device))
     return radius, origin_head, origin_device
 
 
 def _fit_sphere(points, disp='auto'):
-    """Aux function to fit a sphere to an arbitrary set of points"""
+    """Fit a sphere to an arbitrary set of points."""
     from scipy.optimize import fmin_cobyla
     if isinstance(disp, string_types) and disp == 'auto':
         disp = True if logger.level <= 20 else False
@@ -953,7 +987,7 @@ def constraint(center_rad):
 
 
 def _check_origin(origin, info, coord_frame='head', disp=False):
-    """Helper to check or auto-determine the origin"""
+    """Check or auto-determine the origin."""
     if isinstance(origin, string_types):
         if origin != 'auto':
             raise ValueError('origin must be a numerical array, or "auto", '
@@ -1174,7 +1208,7 @@ def read_bem_surfaces(fname, patch_stats=False, s_id=None, verbose=None):
     --------
     write_bem_surfaces, write_bem_solution, make_bem_model
     """
-    from .surface import _complete_surface_info
+    from .surface import complete_surface_info
     # Default coordinate frame
     coord_frame = FIFF.FIFFV_COORD_MRI
     # Open the file, create directory
@@ -1213,7 +1247,7 @@ def read_bem_surfaces(fname, patch_stats=False, s_id=None, verbose=None):
             logger.info('    %d BEM surfaces read' % len(surf))
         if patch_stats:
             for this in surf:
-                _complete_surface_info(this)
+                complete_surface_info(this, copy=False)
     return surf[0] if s_id is not None else surf
 
 

mne/channels/layout.py

@@ -15,7 +15,7 @@
 
 import numpy as np
 
-from ..transforms import _polar_to_cartesian, _cartesian_to_sphere
+from ..transforms import _pol_to_cart, _cart_to_sph
 from ..bem import fit_sphere_to_headshape
 from ..io.pick import pick_types
 from ..io.constants import FIFF
@@ -721,11 +721,8 @@ def _auto_topomap_coords(info, picks, ignore_overlap=False):
         raise ValueError('The following electrodes have overlapping positions:'
                          '\n    ' + str(problematic_electrodes) + '\nThis '
                          'causes problems during visualization.')
-
-    x, y, z = locs3d.T
-    az, el, r = _cartesian_to_sphere(x, y, z)
-    locs2d = np.c_[_polar_to_cartesian(az, np.pi / 2 - el)]
-    return locs2d
+    # use spherical (theta, pol) as (r, theta) for polar->cartesian
+    return _pol_to_cart(_cart_to_sph(locs3d)[:, 1:][:, ::-1])
 
 
 def _topo_to_sphere(pos, eegs):