diff --git a/doc/whats_new.rst b/doc/whats_new.rst index 39afb086670..c1695f0c15d 100644 --- a/doc/whats_new.rst +++ b/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`_ diff --git a/examples/inverse/plot_mixed_norm_inverse.py b/examples/inverse/plot_mixed_norm_inverse.py index 0242fcc05f9..de8cb0af8a4 100644 --- a/examples/inverse/plot_mixed_norm_inverse.py +++ b/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. diff --git a/examples/inverse/plot_time_frequency_mixed_norm_inverse.py b/examples/inverse/plot_time_frequency_mixed_norm_inverse.py index e8739c33753..80be86f0c3e 100644 --- a/examples/inverse/plot_time_frequency_mixed_norm_inverse.py +++ b/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 diff --git a/mne/forward/forward.py b/mne/forward/forward.py index 1cfc3876deb..526421fcfe6 100644 --- a/mne/forward/forward.py +++ b/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. ' diff --git a/mne/inverse_sparse/_gamma_map.py b/mne/inverse_sparse/_gamma_map.py index 567ca6fc99c..a77e386ed54 100644 --- a/mne/inverse_sparse/_gamma_map.py +++ b/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: diff --git a/mne/inverse_sparse/mxne_inverse.py b/mne/inverse_sparse/mxne_inverse.py index 6de1681555f..756e845ec9e 100644 --- a/mne/inverse_sparse/mxne_inverse.py +++ b/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 diff --git a/mne/inverse_sparse/tests/test_gamma_map.py b/mne/inverse_sparse/tests/test_gamma_map.py index 90d225138f7..ff9dee6d519 100644 --- a/mne/inverse_sparse/tests/test_gamma_map.py +++ b/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() diff --git a/mne/inverse_sparse/tests/test_mxne_inverse.py b/mne/inverse_sparse/tests/test_mxne_inverse.py index d239938431f..540c88978a3 100644 --- a/mne/inverse_sparse/tests/test_mxne_inverse.py +++ b/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 diff --git a/tutorials/plot_brainstorm_phantom_elekta.py b/tutorials/plot_brainstorm_phantom_elekta.py index 6806b11920d..663fe7ff3d4 100644 --- a/tutorials/plot_brainstorm_phantom_elekta.py +++ b/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.))