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 <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
