[med-svn] [python-mne] 61/353: ENH : first attempt to make make_inverse_operator work with fixed orientation (small diff with C code to investigate)
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:31 UTC 2015
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 <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list