This is an automated email from the git hooks/post-receive script. yoh pushed a commit to tag 0.4 in repository python-mne.
commit de324e3e8a24817bf1d4a9ec5533a905dea5f47c Author: Alexandre Gramfort <[email protected]> Date: Sun Jan 22 14:54:02 2012 +0100 ENH : first attempt to make make_inverse_operator work with fixed orientation (small diff with C code to investigate) --- mne/forward.py | 11 +++++++++++ mne/minimum_norm/inverse.py | 33 +++++++++++++++++-------------- mne/minimum_norm/tests/test_inverse.py | 36 ++++++++++++++++++++++++++++++++-- 3 files changed, 63 insertions(+), 17 deletions(-) diff --git a/mne/forward.py b/mne/forward.py index 61d017d..ab2520b 100644 --- a/mne/forward.py +++ b/mne/forward.py @@ -432,3 +432,14 @@ def compute_depth_prior(G, exp=0.8, limit=10.0): wpp = (wp / wmax) ** exp depth_prior = np.ravel(wpp[:, None] * np.ones((1, 3))) return depth_prior + + +def compute_depth_prior_fixed(G, exp=0.8, limit=10.0): + """Compute weighting for depth prior for fixed orientation lead field + """ + d = np.sum(G ** 2, axis=0) + w = 1.0 / d + wmax = np.min(w) * (limit ** 2) + wp = np.minimum(w, wmax) + depth_prior = (wp / wmax) ** exp + return depth_prior diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index ece3ad1..101d7b5 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -18,7 +18,7 @@ from ..fiff.tree import dir_tree_find from ..fiff.pick import pick_channels from ..cov import read_cov, prepare_noise_cov -from ..forward import compute_depth_prior +from ..forward import compute_depth_prior, compute_depth_prior_fixed from ..source_space import read_source_spaces_from_tree, \ find_source_space_hemi, _get_vertno from ..transforms import invert_transform, transform_source_space_to @@ -865,12 +865,10 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): # Handle depth prior scaling depth_prior = np.ones(n_dipoles, dtype=gain.dtype) if depth is not None: - if not is_fixed_ori: - depth_prior = compute_depth_prior(gain, exp=depth) + if is_fixed_ori: + depth_prior = compute_depth_prior_fixed(gain, exp=depth) else: - # XXX : how to handle depth_prior with fixed orientation? - warnings.warn('depth_prior is not supported for fixed orientation' - ' forward solutions.') + depth_prior = compute_depth_prior(gain, exp=depth) print "Computing inverse operator with %d channels." % len(ch_names) @@ -878,13 +876,19 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): print 'Whitening lead field matrix.' gain = np.dot(W, gain) - # apply loose orientations - orient_prior = np.ones(n_dipoles, dtype=gain.dtype) - if loose is not None: - print 'Applying loose dipole orientations. Loose value of %s.' % loose - orient_prior[np.mod(np.arange(n_dipoles), 3) != 2] *= loose + source_cov = depth_prior.copy() + depth_prior = dict(data=depth_prior) - source_cov = orient_prior * depth_prior + # apply loose orientations + if not is_fixed_ori: + orient_prior = np.ones(n_dipoles, dtype=gain.dtype) + if loose is not None: + print 'Applying loose dipole orientations. Loose value of %s.' % loose + orient_prior[np.mod(np.arange(n_dipoles), 3) != 2] *= loose + source_cov *= orient_prior + orient_prior = dict(data=orient_prior) + else: + orient_prior = None # Adjusting Source Covariance matrix to make trace of G*R*G' equal # to number of sensors. @@ -896,6 +900,8 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): source_cov *= scaling_source_cov gain *= sqrt(scaling_source_cov) + source_cov = dict(data=source_cov) + # now np.trace(np.dot(gain, gain.T)) == n_nzero # print np.trace(np.dot(gain, gain.T)), n_nzero @@ -904,9 +910,6 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): eigen_fields = dict(data=eigen_fields.T, col_names=ch_names) eigen_leads = dict(data=eigen_leads.T, nrow=eigen_leads.shape[1]) - depth_prior = dict(data=depth_prior) - orient_prior = dict(data=orient_prior) - source_cov = dict(data=source_cov) nave = 1.0 inv_op = dict(eigen_fields=eigen_fields, eigen_leads=eigen_leads, diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 4fdc9c4..567a016 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -2,6 +2,7 @@ import os.path as op import numpy as np from numpy.testing import assert_array_almost_equal, assert_equal from nose.tools import assert_true +import nose from ...datasets import sample from ...label import read_label, label_sign_flip @@ -18,6 +19,8 @@ data_path = sample.data_path(examples_folder) fname_inv = op.join(data_path, 'MEG', 'sample', # 'sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif') 'sample_audvis-meg-oct-6-meg-inv.fif') +fname_inv_fixed = op.join(data_path, 'MEG', 'sample', + 'sample_audvis-meg-oct-6-meg-fixed-inv.fif') fname_vol_inv = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg-vol-7-meg-inv.fif') fname_data = op.join(data_path, 'MEG', 'sample', @@ -35,7 +38,9 @@ label = 'Aud-lh' fname_label = op.join(data_path, 'MEG', 'sample', 'labels', '%s.label' % label) inverse_operator = read_inverse_operator(fname_inv) +inverse_operator_fixed = read_inverse_operator(fname_inv_fixed) label = read_label(fname_label) +noise_cov = Covariance(fname_cov) raw = fiff.Raw(fname_raw) snr = 3.0 lambda2 = 1.0 / snr ** 2 @@ -65,7 +70,6 @@ def test_inverse_operator(): assert_true(stc.data.mean() > 0.1) # Test MNE inverse computation starting from forward operator - noise_cov = Covariance(fname_cov) evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) fwd_op = read_forward_solution(fname_fwd, surf_ori=True) my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov, @@ -77,6 +81,35 @@ def test_inverse_operator(): assert_array_almost_equal(stc.data, my_stc.data, 2) +def test_make_inverse_operator_fixed(): + """Test MNE inverse computation with fixed orientation""" + # XXX : should be fixed and not skipped + raise nose.SkipTest("XFailed Test") + + evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) + fwd_op = read_forward_solution(fname_fwd, force_fixed=True) + inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov, depth=0.8, + loose=None) + + assert_array_almost_equal(inverse_operator_fixed['depth_prior']['data'], + inv_op['depth_prior']['data']) + assert_equal(inverse_operator_fixed['orient_prior'], + inv_op['orient_prior']) + assert_array_almost_equal(inverse_operator_fixed['source_cov']['data'], + inv_op['source_cov']['data']) + + stc_fixed = apply_inverse(evoked, inverse_operator_fixed, lambda2, dSPM) + my_stc = apply_inverse(evoked, inv_op, lambda2, dSPM) + + assert_equal(stc_fixed.times, my_stc.times) + assert_array_almost_equal(stc_fixed.data, my_stc.data, 2) + + # assert_array_almost_equal(inverse_operator_fixed['eigen_fields']['data'], + # inv_op['eigen_fields']['data']) + # assert_array_almost_equal(inverse_operator_fixed['eigen_leads']['data'], + # inv_op['eigen_leads']['data']) + + def test_inverse_operator_volume(): """Test MNE inverse computation on volume source space""" evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) @@ -122,7 +155,6 @@ def test_apply_mne_inverse_fixed_raw(): # create a fixed-orientation inverse operator fwd = read_forward_solution(fname_fwd, force_fixed=True) - noise_cov = Covariance(fname_cov) inv_op = make_inverse_operator(raw.info, fwd, noise_cov, loose=None, depth=0.8) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
