FIX: Fix field mapping

mne/forward/_field_interpolation.py

@@ -261,6 +261,7 @@ def _make_surface_mapping(info, surf, ch_type='meg', trans=None, mode='fast',
         raise ValueError('mode must be "accurate" or "fast", not "%s"' % mode)
 
     # deal with coordinate frames here -- always go to "head" (easiest)
+    orig_surf = surf
     surf = transform_surface_to(deepcopy(surf), 'head', trans)
     n_jobs = check_n_jobs(n_jobs)
     origin = _check_origin(origin, info)
@@ -316,6 +317,8 @@ def _make_surface_mapping(info, surf, ch_type='meg', trans=None, mode='fast',
     logger.info('Field mapping data ready')
 
     fmd['data'] = _compute_mapping_matrix(fmd, info)
+    # bring the original back, whatever coord frame it was in
+    fmd['surf'] = orig_surf
 
     # Remove some unecessary fields
     del fmd['self_dots']
@@ -418,7 +421,6 @@ def make_field_map(evoked, trans='auto', subject=None, subjects_dir=None,
     for this_type, this_surf in zip(types, surfs):
         this_map = _make_surface_mapping(evoked.info, this_surf, this_type,
                                          trans, n_jobs=n_jobs, origin=origin)
-        this_map['surf'] = this_surf  # XXX : a bit weird...
         surf_maps.append(this_map)
 
     return surf_maps

mne/forward/_lead_dots.py

@@ -233,8 +233,10 @@ def _fast_sphere_dot_r0(r, rr1_orig, rr2s, lr1, lr2s, cosmags1, cosmags2s,
     """
     if w1 is None:  # operating on surface, treat independently
         out_shape = (len(rr2s), len(rr1_orig))
+        sum_axis = 1  # operate along second axis only at the end
     else:
         out_shape = (len(rr2s),)
+        sum_axis = None  # operate on flattened array at the end
     out = np.empty(out_shape)
     rr2 = np.concatenate(rr2s)
     lr2 = np.concatenate(lr2s)
@@ -291,7 +293,7 @@ def _fast_sphere_dot_r0(r, rr1_orig, rr2s, lr1, lr2s, cosmags1, cosmags2s,
     if w1 is not None:
         result *= w1[:, np.newaxis]
     for ii, w2 in enumerate(w2s):
-        out[ii] = np.sum(result[:, offset:offset + len(w2)])
+        out[ii] = np.sum(result[:, offset:offset + len(w2)], axis=sum_axis)
         offset += len(w2)
     return out
 

mne/forward/tests/test_field_interpolation.py

@@ -1,8 +1,8 @@
 import numpy as np
 from os import path as op
 from numpy.polynomial import legendre
-from numpy.testing.utils import (assert_allclose, assert_array_equal,
-                                 assert_array_almost_equal)
+from numpy.testing import (assert_allclose, assert_array_equal, assert_equal,
+                           assert_array_almost_equal)
 from nose.tools import assert_raises, assert_true
 
 from mne.forward import _make_surface_mapping, make_field_map
@@ -15,7 +15,7 @@
 from mne.forward._field_interpolation import _setup_dots
 from mne.surface import get_meg_helmet_surf, get_head_surf
 from mne.datasets import testing
-from mne import read_evokeds
+from mne import read_evokeds, pick_types
 from mne.fixes import partial
 from mne.externals.six.moves import zip
 from mne.utils import run_tests_if_main, slow_test
@@ -171,6 +171,28 @@ def test_make_field_map_meg():
                   subjects_dir=subjects_dir, trans=trans_fname)
 
 
+@testing.requires_testing_data
+def test_make_field_map_meeg():
+    """Test making a M/EEG field map onto helmet & head"""
+    evoked = read_evokeds(evoked_fname, baseline=(-0.2, 0.0))[0]
+    picks = pick_types(evoked.info, meg=True, eeg=True)
+    picks = picks[::10]
+    evoked.pick_channels([evoked.ch_names[p] for p in picks])
+    maps = make_field_map(evoked, trans_fname, subject='sample',
+                          subjects_dir=subjects_dir, n_jobs=1)
+    assert_equal(maps[0]['data'].shape, (642, 6))  # EEG->Head
+    assert_equal(maps[1]['data'].shape, (304, 31))  # MEG->Helmet
+    # reasonable ranges
+    for map_ in maps:
+        assert_true(0.5 < map_['data'].max() < 2)
+        assert_true(-2 < map_['data'].min() < -0.5)
+    # calculated from correct looking mapping on 2015/12/26
+    assert_allclose(np.sqrt(np.sum(maps[0]['data'] ** 2)), 16.6088,
+                    atol=1e-3, rtol=1e-3)
+    assert_allclose(np.sqrt(np.sum(maps[1]['data'] ** 2)), 20.1245,
+                    atol=1e-3, rtol=1e-3)
+
+
 def _setup_args(info):
     """Helper to test_as_meg_type_evoked."""
     coils = _create_meg_coils(info['chs'], 'normal', info['dev_head_t'])