MRG: Fix the computation of the orientation in (..)MxNE and gamma-map (…

…#4633)

* fix ori computation

* fix ori for gamma map

* add 3D plot (dipole pos/ori) to the tutorial

* add tests

* Fix ori computation: change method to compute depth weighting

* fix comment

* fix tutorial example

* tests

* add section to tutorial when plotting

* remove code dupes

* typo

* remove warings in gamma map tests

* remove warings in mxne tests

* fix reg parameters to have the same images as before

* check only 1st ori as they are the same

* remove XXX

* update what's new

doc/whats_new.rst

@@ -202,6 +202,9 @@ BUG
 
     - Fix resampling of events along with raw in :func:`mne.io.Raw` to now take into consideration the value of ``first_samp`` by `Chris Bailey`_
 
+    - Fix depth weighting of sparse solvers (:func:`mne.inverse_sparse.mixed_norm`, :func:`mne.inverse_sparse.tf_mixed_norm` and :func:`mne.inverse_sparse.gamma_map`) with free orientation source spaces to improve orientation estimation by `Alex Gramfort`_ and `Yousra Bekhti`_
+
+
 API
 ~~~
     - Add ``skip_by_annotation`` to :meth:`mne.io.Raw.filter` to process data concatenated with e.g. :func:`mne.concatenate_raws` separately. This parameter will default to the old behavior (treating all data as a single block) in 0.15 but will change to ``skip_by_annotation='edge'``, which will separately filter the concatenated chunks separately, in 0.16. This should help prevent potential problems with filter-induced ringing in concatenated files, by `Eric Larson`_

examples/inverse/plot_mixed_norm_inverse.py

@@ -53,7 +53,7 @@
 
 ###############################################################################
 # Run solver
-alpha = 50  # regularization parameter between 0 and 100 (100 is high)
+alpha = 55  # regularization parameter between 0 and 100 (100 is high)
 loose, depth = 0.2, 0.9  # loose orientation & depth weighting
 n_mxne_iter = 10  # if > 1 use L0.5/L2 reweighted mixed norm solver
 # if n_mxne_iter > 1 dSPM weighting can be avoided.

examples/inverse/plot_time_frequency_mixed_norm_inverse.py

@@ -71,7 +71,7 @@
 # Run solver
 
 # alpha_space regularization parameter is between 0 and 100 (100 is high)
-alpha_space = 50.  # spatial regularization parameter
+alpha_space = 30.  # spatial regularization parameter
 # alpha_time parameter promotes temporal smoothness
 # (0 means no temporal regularization)
 alpha_time = 1.  # temporal regularization parameter

mne/forward/forward.py

@@ -618,7 +618,7 @@ def convert_forward_solution(fwd, surf_ori=False, force_fixed=False,
     if use_cps is None:
         if force_fixed:
             use_cps = False
-            warn('The default settings controlling the the application of '
+            warn('The default settings controlling the application of '
                  'cortical patch statistics (cps) in the creation of forward '
                  'operators with fixed orientation will be modified in 0.16. '
                  'The cps (if available) will then be applied by default. '

mne/inverse_sparse/_gamma_map.py

@@ -121,7 +121,8 @@ def denom_fun(x):
             if denom is None:
                 gammas = numer
             else:
-                gammas = numer / denom_fun(denom)
+                gammas = numer / np.maximum(denom_fun(denom),
+                                            np.finfo('float').eps)
         else:
             numer_comb = np.sum(numer.reshape(-1, group_size), axis=1)
             if denom is None:

mne/inverse_sparse/mxne_inverse.py

@@ -6,7 +6,8 @@
 import numpy as np
 from scipy import linalg, signal
 
-from ..source_estimate import SourceEstimate, VolSourceEstimate
+from ..source_estimate import (SourceEstimate, VolSourceEstimate,
+                               _BaseSourceEstimate)
 from ..minimum_norm.inverse import (combine_xyz, _prepare_forward,
                                     _check_reference, _check_loose_forward)
 from ..forward import (compute_orient_prior, is_fixed_orient,
@@ -24,8 +25,7 @@
 @verbose
 def _prepare_weights(forward, gain, source_weighting, weights, weights_min):
     mask = None
-    if isinstance(weights, SourceEstimate):
-        # weights = np.sqrt(np.sum(weights.data ** 2, axis=1))
+    if isinstance(weights, _BaseSourceEstimate):
         weights = np.max(np.abs(weights.data), axis=1)
     weights_max = np.max(weights)
     if weights_min > weights_max:
@@ -61,16 +61,23 @@ def _prepare_gain_column(forward, info, noise_cov, pca, depth, loose, weights,
 
     logger.info('Whitening lead field matrix.')
     gain = np.dot(whitener, gain)
+    is_fixed_ori = is_fixed_orient(forward)
 
     if depth is not None:
-        depth_prior = np.sum(gain ** 2, axis=0) ** depth
+        depth_prior = np.sum(gain ** 2, axis=0)
+        if not is_fixed_ori:
+            depth_prior = depth_prior.reshape(-1, 3).sum(axis=1)
         # Spherical leadfield can be zero at the center
-        depth_prior[depth_prior == 0.] = np.min(depth_prior[depth_prior != 0.])
+        depth_prior[depth_prior == 0.] = np.min(
+            depth_prior[depth_prior != 0.])
+        depth_prior **= depth
+        if not is_fixed_ori:
+            depth_prior = np.repeat(depth_prior, 3)
         source_weighting = np.sqrt(1. / depth_prior)
     else:
         source_weighting = np.ones(gain.shape[1], dtype=gain.dtype)
 
-    assert (is_fixed_orient(forward) or (0 <= loose <= 1))
+    assert (is_fixed_ori or (0 <= loose <= 1))
     if loose is not None and loose < 1.:
         source_weighting *= np.sqrt(compute_orient_prior(forward, loose))
 
@@ -178,7 +185,6 @@ def _make_sparse_stc(X, active_set, forward, tmin, tstep,
 @verbose
 def _make_dipoles_sparse(X, active_set, forward, tmin, tstep, M, M_est,
                          active_is_idx=False, verbose=None):
-
     times = tmin + tstep * np.arange(X.shape[1])
 
     if not active_is_idx:
@@ -209,10 +215,12 @@ def _make_dipoles_sparse(X, active_set, forward, tmin, tstep, M, M_est,
             if forward['surf_ori']:
                 X_ = np.dot(forward['source_nn'][i_dip *
                             n_dip_per_pos:(i_dip + 1) * n_dip_per_pos].T, X_)
+
             amplitude = np.sqrt(np.sum(X_ ** 2, axis=0))
             i_ori = np.zeros((len(times), 3))
             i_ori[amplitude > 0.] = (X_[:, amplitude > 0.] /
                                      amplitude[amplitude > 0.]).T
+
         dipoles.append(Dipole(times, i_pos, amplitude, i_ori, gof))
 
     return dipoles

mne/inverse_sparse/tests/test_gamma_map.py

@@ -124,5 +124,22 @@ def test_gamma_map_vol_sphere():
 
     assert_array_almost_equal(stc.times, evoked.times, 5)
 
+    # Compare orientation obtained using fit_dipole and gamma_map
+    # for a simulated evoked containing a single dipole
+    stc = mne.VolSourceEstimate(50e-9 * np.random.RandomState(42).randn(1, 4),
+                                vertices=stc.vertices[:1],
+                                tmin=stc.tmin,
+                                tstep=stc.tstep)
+    evoked_dip = mne.simulation.simulate_evoked(fwd, stc, info, cov, nave=1e9,
+                                                use_cps=True)
+
+    dip_gmap = gamma_map(evoked_dip, fwd, cov, 0.1, return_as_dipoles=True)
+
+    amp_max = [np.max(d.amplitude) for d in dip_gmap]
+    dip_gmap = dip_gmap[np.argmax(amp_max)]
+    assert_true(dip_gmap[0].pos[0] in src[0]['rr'][stc.vertices])
+
+    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
+    assert_true(np.abs(np.dot(dip_fit.ori[0], dip_gmap.ori[0])) > 0.99)
 
 run_tests_if_main()

mne/inverse_sparse/tests/test_mxne_inverse.py

@@ -170,6 +170,26 @@ def test_mxne_vol_sphere():
     assert_true(isinstance(stc, VolSourceEstimate))
     assert_array_almost_equal(stc.times, evoked_l21.times, 5)
 
+    # Compare orientation obtained using fit_dipole and gamma_map
+    # for a simulated evoked containing a single dipole
+    stc = mne.VolSourceEstimate(50e-9 * np.random.RandomState(42).randn(1, 4),
+                                vertices=stc.vertices[:1],
+                                tmin=stc.tmin,
+                                tstep=stc.tstep)
+    evoked_dip = mne.simulation.simulate_evoked(fwd, stc, info, cov, nave=1e9,
+                                                use_cps=True)
+
+    dip_mxne = mixed_norm(evoked_dip, fwd, cov, alpha=80,
+                          n_mxne_iter=1, maxit=30, tol=1e-8,
+                          active_set_size=10, return_as_dipoles=True)
+
+    amp_max = [np.max(d.amplitude) for d in dip_mxne]
+    dip_mxne = dip_mxne[np.argmax(amp_max)]
+    assert_true(dip_mxne.pos[0] in src[0]['rr'][stc.vertices])
+
+    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
+    assert_true(np.abs(np.dot(dip_fit.ori[0], dip_mxne.ori[0])) > 0.99)
+
     # Do with TF-MxNE for test memory savings
     alpha_space = 60.  # spatial regularization parameter
     alpha_time = 1.  # temporal regularization parameter

tutorials/plot_brainstorm_phantom_elekta.py

@@ -30,6 +30,7 @@
 from mne.datasets.brainstorm import bst_phantom_elekta
 from mne.io import read_raw_fif
 
+from mayavi import mlab
 print(__doc__)
 
 ###############################################################################
@@ -137,3 +138,22 @@
 
 fig.tight_layout()
 plt.show()
+
+
+###############################################################################
+# Let's plot the positions and the orientations of the actual and the estimated
+# dipoles
+def plot_pos_ori(pos, ori, color=(0., 0., 0.)):
+    mlab.points3d(pos[:, 0], pos[:, 1], pos[:, 2], scale_factor=0.005,
+                  color=color)
+    mlab.quiver3d(pos[:, 0], pos[:, 1], pos[:, 2],
+                  ori[:, 0], ori[:, 1], ori[:, 2],
+                  scale_factor=0.03,
+                  color=color)
+
+mne.viz.plot_alignment(evoked.info, bem=sphere, surfaces=[])
+
+# Plot the position and the orientation of the actual dipole
+plot_pos_ori(actual_pos, actual_ori, color=(1., 0., 0.))
+# Plot the position and the orientation of the estimated dipole
+plot_pos_ori(dip.pos, dip.ori, color=(0., 0., 1.))