[med-svn] [python-mne] 61/376: API : do not read bad channels from forward operator ENH : adding possibility to whiten forward solutions using noise covariance matrix.
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:07 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to annotated tag v0.1
in repository python-mne.
commit 4a38978ba3aab0884c690db31f18c547ee757d70
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Jan 31 17:14:05 2011 -0500
API : do not read bad channels from forward operator
ENH : adding possibility to whiten forward solutions using noise covariance
matrix.
---
examples/{read_forward.py => plot_read_forward.py} | 14 ++--
mne/cov.py | 49 +++++++++++++-
mne/fiff/pick.py | 61 +++++++++++++++--
mne/forward.py | 77 ++++++++--------------
4 files changed, 139 insertions(+), 62 deletions(-)
diff --git a/examples/read_forward.py b/examples/plot_read_forward.py
similarity index 77%
rename from examples/read_forward.py
rename to examples/plot_read_forward.py
index 583f6d3..caa8338 100644
--- a/examples/read_forward.py
+++ b/examples/plot_read_forward.py
@@ -15,8 +15,8 @@ import mne
fname = os.environ['MNE_SAMPLE_DATASET_PATH']
fname += '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
-data = mne.read_forward_solution(fname)
-leadfield = data['sol']['data']
+fwd = mne.read_forward_solution(fname)
+leadfield = fwd['sol']['data']
print "Leadfield size : %d x %d" % leadfield.shape
@@ -24,17 +24,17 @@ print "Leadfield size : %d x %d" % leadfield.shape
# Show result
import pylab as pl
-pl.matshow(leadfield[:306,:500])
+pl.matshow(leadfield[:,:500])
pl.xlabel('sources')
pl.ylabel('sensors')
pl.title('Lead field matrix')
pl.show()
# 3D source space
-lh_points = data['src'][0]['rr']
-lh_faces = data['src'][0]['use_tris']
-rh_points = data['src'][1]['rr']
-rh_faces = data['src'][1]['use_tris']
+lh_points = fwd['src'][0]['rr']
+lh_faces = fwd['src'][0]['use_tris']
+rh_points = fwd['src'][1]['rr']
+rh_faces = fwd['src'][1]['use_tris']
from enthought.mayavi import mlab
mlab.triangular_mesh(lh_points[:,0], lh_points[:,1], lh_points[:,2], lh_faces)
mlab.triangular_mesh(rh_points[:,0], rh_points[:,1], rh_points[:,2], rh_faces)
diff --git a/mne/cov.py b/mne/cov.py
index 738e9d1..1aa8288 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -18,7 +18,7 @@ from .fiff.write import start_block, end_block, write_int, write_name_list, \
write_double, write_float_matrix, start_file, end_file
from .fiff.proj import write_proj
from .fiff import fiff_open
-from .fiff.pick import pick_types
+from .fiff.pick import pick_types, pick_channels_forward
class Covariance(object):
@@ -152,6 +152,53 @@ class Covariance(object):
return ave_whiten, W
+ def whiten_evoked_and_forward(self, ave, fwd, eps=0.2):
+ """Whiten an evoked data set and a forward solution
+
+ The whitening matrix is estimated and then multiplied to
+ forward solution a.k.a. the leadfield matrix.
+ It makes the additive white noise assumption of MNE
+ realistic.
+
+ Parameters
+ ----------
+ ave : evoked data
+ A evoked data set read with fiff.read_evoked
+ fwd : forward data
+ A forward solution read with mne.read_forward
+ eps : float
+ The regularization factor used.
+
+ Returns
+ -------
+ ave : evoked data
+ A evoked data set read with fiff.read_evoked
+ fwd : evoked data
+ Forward solution after whitening.
+ W : array of shape [n_channels, n_channels]
+ The whitening matrix
+ """
+ # handle evoked
+ ave_whiten, W = self.whiten_evoked(ave, eps=eps)
+
+ ave_ch_names = [ch['ch_name'] for ch in ave_whiten['info']['chs']]
+
+ # handle forward (keep channels in covariance matrix)
+ fwd_whiten = copy.copy(fwd)
+ ind = [fwd_whiten['sol']['row_names'].index(name)
+ for name in self._cov['names']]
+ fwd_whiten['sol']['data'][ind] = np.dot(W,
+ fwd_whiten['sol']['data'][ind])
+ fwd_whiten['sol']['row_names'] = [fwd_whiten['sol']['row_names'][k]
+ for k in ind]
+ fwd_whiten['chs'] = [fwd_whiten['chs'][k] for k in ind]
+
+ # keep in forward the channels in the evoked dataset
+ fwd_whiten = pick_channels_forward(fwd, include=ave_ch_names,
+ exclude=ave['info']['bads'])
+
+ return ave_whiten, fwd_whiten, W
+
def __repr__(self):
s = "kind : %s" % self.kind
s += ", size : %s x %s" % self.data.shape
diff --git a/mne/fiff/pick.py b/mne/fiff/pick.py
index 08a672a..9bf3e72 100644
--- a/mne/fiff/pick.py
+++ b/mne/fiff/pick.py
@@ -32,7 +32,7 @@ def pick_channels(ch_names, include, exclude):
"""
sel = []
for k, name in enumerate(ch_names):
- if (include is [] or name in include) and name not in exclude:
+ if (len(include) == 0 or name in include) and name not in exclude:
sel.append(k)
sel = np.unique(sel)
np.sort(sel)
@@ -40,7 +40,7 @@ def pick_channels(ch_names, include, exclude):
def pick_types(info, meg=True, eeg=False, stim=False, include=[], exclude=[]):
- """Pick channels
+ """Pick channels by type and names
Parameters
----------
@@ -146,7 +146,7 @@ def pick_channels_evoked(orig, include=[], exclude=[]):
exclude are None it returns orig without copy.
"""
- if include is None and exclude is None:
+ if len(include) == 0 and len(exclude) == 0:
return orig
sel = pick_channels(orig['info']['ch_names'], include=include,
@@ -155,7 +155,7 @@ def pick_channels_evoked(orig, include=[], exclude=[]):
if len(sel) == 0:
raise ValueError, 'Warning : No channels match the selection.'
- res = orig.copy()
+ res = copy(orig)
#
# Modify the measurement info
#
@@ -166,3 +166,56 @@ def pick_channels_evoked(orig, include=[], exclude=[]):
res['evoked']['epochs'] = res['evoked']['epochs'][sel,:]
return res
+
+
+def pick_channels_forward(orig, include=[], exclude=[]):
+ """Pick channels from forward operator
+
+ Parameters
+ ----------
+ orig : dict
+ A forward solution
+
+ include : list of string, (optional)
+ List of channels to include. (if None, include all available)
+
+ exclude : list of string, (optional)
+ Channels to exclude (if None, do not exclude any)
+
+ Returns
+ -------
+ res : dict
+ Evoked data restricted to selected channels. If include and
+ exclude are None it returns orig without copy.
+ """
+
+ if len(include) == 0 and len(exclude) == 0:
+ return orig
+
+ sel = pick_channels(orig['sol']['row_names'], include=include,
+ exclude=exclude)
+
+ fwd = copy(orig)
+
+ # Do we have something?
+ nuse = len(sel)
+ if nuse == 0:
+ raise ValueError, 'Nothing remains after picking'
+
+ print '\t%d out of %d channels remain after picking' % (nuse,
+ fwd['nchan'])
+
+ # Pick the correct rows of the forward operator
+ fwd['nchan'] = nuse
+ fwd['sol']['data'] = fwd['sol']['data'][sel, :]
+ fwd['sol']['nrow'] = nuse
+ fwd['sol']['row_names'] = [fwd['sol']['row_names'][k] for k in sel]
+ fwd['chs'] = [fwd['chs'][k] for k in sel]
+
+ if fwd['sol_grad'] is not None:
+ fwd['sol_grad']['data'] = fwd['sol_grad']['data'][sel, :]
+ fwd['sol_grad']['nrow'] = nuse
+ fwd['sol_grad']['row_names'] = [fwd['sol_grad']['row_names'][k]
+ for k in sel]
+
+ return fwd
diff --git a/mne/forward.py b/mne/forward.py
index f62df1c..1e7c28a 100644
--- a/mne/forward.py
+++ b/mne/forward.py
@@ -9,13 +9,14 @@ from scipy import linalg
from .fiff.constants import FIFF
from .fiff.open import fiff_open
from .fiff.tree import dir_tree_find
-from .fiff.channels import _read_bad_channels
-from .fiff.tag import find_tag
+from .fiff.tag import find_tag, read_tag
from .fiff.matrix import _read_named_matrix, _transpose_named_matrix
+from .fiff.pick import pick_channels_forward
from .source_space import read_source_spaces, find_source_space_hemi
from .transforms import transform_source_space_to, invert_transform
+
def _block_diag(A, n):
"""Constructs a block diagonal from a packed structure
@@ -146,7 +147,7 @@ def _read_one(fid, node):
def read_forward_solution(fname, force_fixed=False, surf_ori=False,
- include=None, exclude=None):
+ include=[], exclude=[]):
"""Read a forward solution a.k.a. lead field
Parameters
@@ -161,10 +162,12 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False,
Use surface based source coordinate system?
include: list, optional
- List of names of channels to include
+ List of names of channels to include. If empty all channels
+ are included.
exclude: list, optional
- List of names of channels to exclude
+ List of names of channels to exclude. If empty include all
+ channels.
Returns
-------
@@ -190,6 +193,21 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False,
raise ValueError, 'No parent MRI information in %s' % fname
parent_mri = parent_mri[0]
+ # Parent MEG data
+ parent_meg = dir_tree_find(tree, FIFF.FIFFB_MNE_PARENT_MEAS_FILE)
+ if len(parent_meg) == 0:
+ fid.close()
+ raise ValueError, 'No parent MEG information in %s' % fname
+ parent_meg = parent_meg[0]
+
+ chs = list()
+ for k in range(parent_meg['nent']):
+ kind = parent_meg['directory'][k].kind
+ pos = parent_meg['directory'][k].pos
+ if kind == FIFF.FIFF_CH_INFO:
+ tag = read_tag(fid, pos)
+ chs.append(tag.data)
+
try:
src = read_source_spaces(fid, False, tree)
except Exception as inst:
@@ -201,11 +219,6 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False,
fwd = None
- # Bad channel list
- bads = _read_bad_channels(fid, tree)
-
- print '\t%d bad channels read' % len(bads)
-
# Locate and read the forward solutions
megnode = None
eegnode = None
@@ -383,46 +396,10 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False,
fwd['source_nn'] = np.kron(np.ones((fwd['nsource'], 1)), np.eye(3))
print '[done]'
- # Do the channel selection
- if include is not None or exclude is not None or bads is not None:
- # First do the channels to be included
- pick = np.ones(fwd['nchan'], dtype=np.bool)
- if include is not None:
- for p in range(len(fwd['sol']['row_names'])):
- if fwd['sol']['row_names'][p] in include:
- pick[p] = True
-
- # Then exclude what needs to be excluded
- if exclude is not None:
- for p in range(len(fwd['sol']['row_names'])):
- if fwd['sol']['row_names'][p] in exclude:
- pick[p] = False
-
- if bads is not None:
- for k in range(len(bads)):
- for p in range(len(fwd['sol']['row_names'])):
- if fwd['sol']['row_names'][p] in bads:
- pick[p] = False
-
- # Do we have something?
- nuse = pick.sum()
- if nuse == 0:
- raise ValueError, 'Nothing remains after picking'
-
- print '\t%d out of %d channels remain after picking' % (nuse,
- fwd['nchan'])
-
- # Pick the correct rows of the forward operator
- fwd['nchan'] = nuse
- fwd['sol']['data'] = fwd['sol']['data'][pick == 1, :]
- fwd['sol']['nrow'] = nuse
- fwd['sol']['row_names'] = [fwd['sol']['row_names'][k]
- for k in range(len(pick)) if pick[k]]
+ # Add channel information
+ fwd['chs'] = chs
- if fwd['sol_grad'] is not None:
- fwd['sol_grad']['data'] = fwd['sol_grad']['data'][pick == 1, :]
- fwd['sol_grad']['nrow'] = nuse
- fwd['sol_grad']['row_names'] = [fwd['sol_grad']['row_names'][k]
- for k in range(len(pick)) if pick[k]]
+ fwd = pick_channels_forward(fwd, include=include, exclude=exclude)
return fwd
+
--
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