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 4f67f889819ef1b593dc178d389c7d5d3122710b Author: Alexandre Gramfort <[email protected]> Date: Thu Jan 27 21:54:26 2011 -0500 ENH : pick by channel type ENH : data whitening with noise covariance matrix ENH : moving transform functions in transforms.py --- examples/plot_read_epochs.py | 38 +++++++---- examples/plot_whitened_evoked_data.py | 56 +++++++++++++++ mne/cov.py | 73 ++++++++++++++++++++ mne/fiff/constants.py | 67 ++++++++++++++++++ mne/fiff/pick.py | 18 +++-- mne/forward.py | 50 +------------- mne/inverse.py | 8 +-- mne/transforms.py | 125 ++++++++++++++++++++++++++++++++++ 8 files changed, 367 insertions(+), 68 deletions(-) diff --git a/examples/plot_read_epochs.py b/examples/plot_read_epochs.py index f1745f4..98cfcbb 100644 --- a/examples/plot_read_epochs.py +++ b/examples/plot_read_epochs.py @@ -41,37 +41,51 @@ events = mne.read_events(event_fname) include = [] # or stim channel ['STI 014'] exclude = raw['info']['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more -# MEG -meg_picks = fiff.pick_types(raw['info'], meg=True, eeg=False, stim=False, +# MEG Magnetometers +meg_mag_picks = fiff.pick_types(raw['info'], meg='mag', eeg=False, stim=False, include=include, exclude=exclude) -meg_data, times, channel_names = mne.read_epochs(raw, events, event_id, - tmin, tmax, picks=meg_picks, baseline=(None, 0)) -meg_epochs = np.array([d['epoch'] for d in meg_data]) # build 3D matrix -meg_evoked_data = np.mean(meg_epochs, axis=0) # compute evoked fields +meg_mag_data, times, channel_names = mne.read_epochs(raw, events, event_id, + tmin, tmax, picks=meg_mag_picks, baseline=(None, 0)) +meg_mag_epochs = np.array([d['epoch'] for d in meg_mag_data]) # as 3D matrix +meg_mag_evoked_data = np.mean(meg_mag_epochs, axis=0) # compute evoked fields + +# MEG +meg_grad_picks = fiff.pick_types(raw['info'], meg='grad', eeg=False, + stim=False, include=include, exclude=exclude) +meg_grad_data, times, channel_names = mne.read_epochs(raw, events, event_id, + tmin, tmax, picks=meg_grad_picks, baseline=(None, 0)) +meg_grad_epochs = np.array([d['epoch'] for d in meg_grad_data]) # as 3D matrix +meg_grad_evoked_data = np.mean(meg_grad_epochs, axis=0) # compute evoked fields # EEG eeg_picks = fiff.pick_types(raw['info'], meg=False, eeg=True, stim=False, include=include, exclude=exclude) eeg_data, times, channel_names = mne.read_epochs(raw, events, event_id, tmin, tmax, picks=eeg_picks, baseline=(None, 0)) -eeg_epochs = np.array([d['epoch'] for d in eeg_data]) # build 3D matrix +eeg_epochs = np.array([d['epoch'] for d in eeg_data]) # as 3D matrix eeg_evoked_data = np.mean(eeg_epochs, axis=0) # compute evoked potentials ############################################################################### # View evoked response import pylab as pl pl.clf() -pl.subplot(2, 1, 1) -pl.plot(times, meg_evoked_data.T) +pl.subplot(3, 1, 1) +pl.plot(times, meg_mag_evoked_data.T) pl.xlim([times[0], times[-1]]) pl.xlabel('time (ms)') pl.ylabel('Magnetic Field (T)') -pl.title('MEG evoked field') -pl.subplot(2, 1, 2) +pl.title('MEG (Magnetometers) evoked field') +pl.subplot(3, 1, 2) +pl.plot(times, meg_grad_evoked_data.T) +pl.xlim([times[0], times[-1]]) +pl.xlabel('time (ms)') +pl.ylabel('Magnetic Field (T/m)') +pl.title('MEG (Gradiometers) evoked field') +pl.subplot(3, 1, 3) pl.plot(times, eeg_evoked_data.T) pl.xlim([times[0], times[-1]]) pl.xlabel('time (ms)') pl.ylabel('Potential (V)') pl.title('EEG evoked potential') -pl.subplots_adjust(0.175, 0.04, 0.94, 0.94, 0.2, 0.33) +pl.subplots_adjust(0.175, 0.04, 0.94, 0.94, 0.2, 0.53) pl.show() diff --git a/examples/plot_whitened_evoked_data.py b/examples/plot_whitened_evoked_data.py new file mode 100644 index 0000000..d863857 --- /dev/null +++ b/examples/plot_whitened_evoked_data.py @@ -0,0 +1,56 @@ +""" +===================================================== +Whiten an evoked date using a noise covariance matrix +===================================================== + +""" +# Author: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import os +import mne +from mne import fiff + +fname = os.environ['MNE_SAMPLE_DATASET_PATH'] +fname += '/MEG/sample/sample_audvis-ave.fif' +cov_fname = os.environ['MNE_SAMPLE_DATASET_PATH'] +cov_fname += '/MEG/sample/sample_audvis-cov.fif' + +# Reading +ave = fiff.read_evoked(fname, setno=0, baseline=(None, 0)) +cov = mne.Covariance() +cov.load(cov_fname) + +ave_whiten, W = cov.whiten_evoked(ave) + +bads = ave_whiten['info']['bads'] +ind_meg_grad = fiff.pick_types(ave['info'], meg='grad', exclude=bads) +ind_meg_mag = fiff.pick_types(ave['info'], meg='mag', exclude=bads) +ind_eeg = fiff.pick_types(ave['info'], meg=False, eeg=True, exclude=bads) + +############################################################################### +# Show result +import pylab as pl +pl.clf() +pl.subplot(3, 1, 1) +pl.plot(ave['evoked']['times'], + ave_whiten['evoked']['epochs'][ind_meg_grad,:].T) +pl.title('MEG Planar Gradiometers') +pl.xlabel('time (s)') +pl.ylabel('MEG data') +pl.subplot(3, 1, 2) +pl.plot(ave['evoked']['times'], + ave_whiten['evoked']['epochs'][ind_meg_mag,:].T) +pl.title('MEG Magnetometers') +pl.xlabel('time (s)') +pl.ylabel('MEG data') +pl.subplot(3, 1, 3) +pl.plot(ave['evoked']['times'], ave_whiten['evoked']['epochs'][ind_eeg,:].T) +pl.title('EEG') +pl.xlabel('time (s)') +pl.ylabel('EEG data') +pl.subplots_adjust(0.1, 0.08, 0.94, 0.94, 0.2, 0.63) +pl.show() diff --git a/mne/cov.py b/mne/cov.py index 4439aae..738e9d1 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -4,7 +4,9 @@ # License: BSD (3-clause) import os +import copy import numpy as np +from scipy import linalg from .fiff.constants import FIFF from .fiff.tag import find_tag @@ -16,6 +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 class Covariance(object): @@ -79,6 +82,76 @@ class Covariance(object): self.data = cov / n_samples # XXX : check print '[done]' + def _regularize(self, data, variances, ch_names, eps): + """Operates inplace in data + """ + if len(ch_names) > 0: + ind = [self._cov['names'].index(name) for name in ch_names] + reg = eps * np.mean(variances[ind]) + for ii in ind: + data[ind,ind] += reg + + def whiten_evoked(self, ave, eps=0.2): + """Whiten an evoked data file + + The whitening matrix is estimated and then multiplied to data. + It makes the additive white noise assumption of MNE + realistic. + + Parameters + ---------- + ave : evoked data + A evoked data set read with fiff.read_evoked + eps : float + The regularization factor used. + + Returns + ------- + ave : evoked data + Evoked data set after whitening. + W : array of shape [n_channels, n_channels] + The whitening matrix + """ + + data = self.data.copy() # will be the regularized covariance + variances = np.diag(data) + + # Add (eps x identity matrix) to magnetometers only. + # This is based on the mean magnetometer variance like MNE C-code does it. + mag_ind = pick_types(ave['info'], meg='mag', eeg=False, stim=False) + mag_names = [ave['info']['chs'][k]['ch_name'] for k in mag_ind] + self._regularize(data, variances, mag_names, eps) + + # Add (eps x identity matrix) to gradiometers only. + grad_ind = pick_types(ave['info'], meg='grad', eeg=False, stim=False) + grad_names = [ave['info']['chs'][k]['ch_name'] for k in grad_ind] + self._regularize(data, variances, grad_names, eps) + + # Add (eps x identity matrix) to eeg only. + eeg_ind = pick_types(ave['info'], meg=False, eeg=True, stim=False) + eeg_names = [ave['info']['chs'][k]['ch_name'] for k in eeg_ind] + self._regularize(data, variances, eeg_names, eps) + + d, V = linalg.eigh(data) # Compute eigen value decomposition. + + # Compute the unique square root inverse, which is a whitening matrix. + # This matrix can be multiplied with data and leadfield matrix to get + # whitened inverse solutions. + d = 1.0 / np.sqrt(d) + W = np.dot(V, d[:,None] * V.T) + + # Get all channel indices + n_channels = len(ave['info']['chs']) + ave_ch_names = [ave['info']['chs'][k]['ch_name'] + for k in range(n_channels)] + ind = [ave_ch_names.index(name) for name in self._cov['names']] + + ave_whiten = copy.copy(ave) + ave_whiten['evoked']['epochs'][ind] = np.dot(W, + ave['evoked']['epochs'][ind]) + + return ave_whiten, W + def __repr__(self): s = "kind : %s" % self.kind s += ", size : %s x %s" % self.data.shape diff --git a/mne/fiff/constants.py b/mne/fiff/constants.py index b312aa5..2af0acf 100644 --- a/mne/fiff/constants.py +++ b/mne/fiff/constants.py @@ -369,3 +369,70 @@ FIFF.FIFFT_CH_POS_STRUCT = 34 FIFF.FIFFT_COORD_TRANS_STRUCT = 35 FIFF.FIFFT_DIG_STRING_STRUCT = 36 FIFF.FIFFT_STREAM_SEGMENT_STRUCT = 37 +# +# Units of measurement +# +FIFF.FIFF_UNIT_NONE = -1 +# +# SI base units +# +FIFF.FIFF_UNIT_M = 1 +FIFF.FIFF_UNIT_KG = 2 +FIFF.FIFF_UNIT_SEC = 3 +FIFF.FIFF_UNIT_A = 4 +FIFF.FIFF_UNIT_K = 5 +FIFF.FIFF_UNIT_MOL = 6 +# +# SI Supplementary units +# +FIFF.FIFF_UNIT_RAD = 7 +FIFF.FIFF_UNIT_SR = 8 +# +# SI base candela +# +FIFF.FIFF_UNIT_CD = 9 +# +# SI derived units +# +FIFF.FIFF_UNIT_HZ = 101 +FIFF.FIFF_UNIT_N = 102 +FIFF.FIFF_UNIT_PA = 103 +FIFF.FIFF_UNIT_J = 104 +FIFF.FIFF_UNIT_W = 105 +FIFF.FIFF_UNIT_C = 106 +FIFF.FIFF_UNIT_V = 107 +FIFF.FIFF_UNIT_F = 108 +FIFF.FIFF_UNIT_OHM = 109 +FIFF.FIFF_UNIT_MHO = 110 +FIFF.FIFF_UNIT_WB = 111 +FIFF.FIFF_UNIT_T = 112 +FIFF.FIFF_UNIT_H = 113 +FIFF.FIFF_UNIT_CEL = 114 +FIFF.FIFF_UNIT_LM = 115 +FIFF.FIFF_UNIT_LX = 116 +# +# Others we need +# +FIFF.FIFF_UNIT_T_M = 201 # T/m +FIFF.FIFF_UNIT_AM = 202 # Am +FIFF.FIFF_UNIT_AM_M2 = 203 # Am/m^2 +FIFF.FIFF_UNIT_AM_M3 = 204 # Am/m^3 +# +# Multipliers +# +FIFF.FIFF_UNITM_E = 18 +FIFF.FIFF_UNITM_PET = 15 +FIFF.FIFF_UNITM_T = 12 +FIFF.FIFF_UNITM_MEG = 6 +FIFF.FIFF_UNITM_K = 3 +FIFF.FIFF_UNITM_H = 2 +FIFF.FIFF_UNITM_DA = 1 +FIFF.FIFF_UNITM_NONE = 0 +FIFF.FIFF_UNITM_D = -1 +FIFF.FIFF_UNITM_C = -2 +FIFF.FIFF_UNITM_M = -3 +FIFF.FIFF_UNITM_MU = -6 +FIFF.FIFF_UNITM_N = -9 +FIFF.FIFF_UNITM_P = -12 +FIFF.FIFF_UNITM_F = -15 +FIFF.FIFF_UNITM_A = -18 diff --git a/mne/fiff/pick.py b/mne/fiff/pick.py index 2bbc266..08a672a 100644 --- a/mne/fiff/pick.py +++ b/mne/fiff/pick.py @@ -47,8 +47,10 @@ def pick_types(info, meg=True, eeg=False, stim=False, include=[], exclude=[]): info : dict The measurement info - meg : bool - Is True include MEG channels + meg : bool or string + Is True include MEG channels or False include None + If string it can be 'mag' or 'grad' to select only gradiometers + or magnetometers. eeg : bool Is True include EEG channels @@ -72,9 +74,15 @@ def pick_types(info, meg=True, eeg=False, stim=False, include=[], exclude=[]): for k in range(nchan): kind = info['chs'][k].kind - if (kind == FIFF.FIFFV_MEG_CH or kind == FIFF.FIFFV_REF_MEG_CH) \ - and meg: - pick[k] = True + if (kind == FIFF.FIFFV_MEG_CH or kind == FIFF.FIFFV_REF_MEG_CH): + if meg == True: + pick[k] = True + elif (meg is 'grad' + and info['chs'][k]['unit'] == FIFF.FIFF_UNIT_T_M): + pick[k] = True + elif (meg is 'mag' + and info['chs'][k]['unit'] == FIFF.FIFF_UNIT_T): + pick[k] = True elif kind == FIFF.FIFFV_EEG_CH and eeg: pick[k] = True elif kind == FIFF.FIFFV_STIM_CH and stim: diff --git a/mne/forward.py b/mne/forward.py index ed5acee..f62df1c 100644 --- a/mne/forward.py +++ b/mne/forward.py @@ -3,7 +3,6 @@ # # License: BSD (3-clause) -import copy import numpy as np from scipy import linalg @@ -15,7 +14,7 @@ from .fiff.tag import find_tag from .fiff.matrix import _read_named_matrix, _transpose_named_matrix 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 @@ -83,49 +82,6 @@ def _block_diag(A, n): return bd -def _transform_source_space_to(src, dest, trans): - """ - % - % [res] = mne_transform_source_space_to(src,dest,trans) - % - % Transform source space data to the desired coordinate system - % - % fname - The name of the file - % include - Include these channels (optional) - % exclude - Exclude these channels (optional) - % - """ - - if src['coord_frame'] == dest: - res = src - return res - - if trans['to'] == src['coord_frame'] and trans['from_'] == dest: - trans = _invert_transform(trans) - elif trans['from_'] != src['coord_frame'] or trans['to'] != dest: - raise ValueError, 'Cannot transform the source space using this ' \ - 'coordinate transformation' - - t = trans['trans'][:3,:] - res = src - res['coord_frame'] = dest - - res['rr'] = np.dot(np.c_[res['rr'], np.ones((res['np'], 1))], t.T) - res['nn'] = np.dot(np.c_[res['nn'], np.zeros((res['np'], 1))], t.T) - return res - - -def _invert_transform(trans): - """Invert a transformation between coordinate systems - """ - itrans = copy.copy(trans) - aux = itrans['from_'] - itrans['from_'] = itrans['to'] - itrans['to'] = aux - itrans['trans'] = linalg.inv(itrans['trans']) - return itrans - - def _read_one(fid, node): """Read all interesting stuff for one forward solution """ @@ -327,7 +283,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, mri_head_t = tag.data; if (mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): - mri_head_t = _invert_transform(mri_head_t) + mri_head_t = invert_transform(mri_head_t) if (mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): fid.close() @@ -349,7 +305,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, nuse = 0 for s in src: try: - s = _transform_source_space_to(s, fwd['coord_frame'], mri_head_t) + s = transform_source_space_to(s, fwd['coord_frame'], mri_head_t) except Exception as inst: raise ValueError, 'Could not transform source space (%s)' % inst diff --git a/mne/inverse.py b/mne/inverse.py index 3a22873..c1d1566 100644 --- a/mne/inverse.py +++ b/mne/inverse.py @@ -17,8 +17,8 @@ from .fiff.pick import pick_channels_evoked from .cov import read_cov from .source_space import read_source_spaces, find_source_space_hemi -from .forward import _invert_transform, _transform_source_space_to, \ - _block_diag +from .forward import _block_diag +from .transforms import invert_transform, transform_source_space_to def read_inverse_operator(fname): @@ -199,7 +199,7 @@ def read_inverse_operator(fname): mri_head_t = tag.data if mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or \ mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD: - mri_head_t = _invert_transform(mri_head_t) + mri_head_t = invert_transform(mri_head_t) if mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or \ mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD: fid.close() @@ -237,7 +237,7 @@ def read_inverse_operator(fname): nuse = 0 for k in range(len(inv['src'])): try: - inv['src'][k] = _transform_source_space_to(inv['src'][k], + inv['src'][k] = transform_source_space_to(inv['src'][k], inv['coord_frame'], mri_head_t) except Exception as inst: fid.close() diff --git a/mne/transforms.py b/mne/transforms.py new file mode 100644 index 0000000..ba6323c --- /dev/null +++ b/mne/transforms.py @@ -0,0 +1,125 @@ +import copy +import numpy as np +from scipy import linalg + +from .fiff import FIFF + + +def invert_transform(trans): + """Invert a transformation between coordinate systems + """ + itrans = copy.copy(trans) + aux = itrans['from_'] + itrans['from_'] = itrans['to'] + itrans['to'] = aux + itrans['trans'] = linalg.inv(itrans['trans']) + return itrans + + +def transform_source_space_to(src, dest, trans): + """ + % + % [res] = mne_transform_source_space_to(src,dest,trans) + % + % Transform source space data to the desired coordinate system + % + % fname - The name of the file + % include - Include these channels (optional) + % exclude - Exclude these channels (optional) + % + + XXX + """ + + if src['coord_frame'] == dest: + res = src + return res + + if trans['to'] == src['coord_frame'] and trans['from_'] == dest: + trans = invert_transform(trans) + elif trans['from_'] != src['coord_frame'] or trans['to'] != dest: + raise ValueError, 'Cannot transform the source space using this ' \ + 'coordinate transformation' + + t = trans['trans'][:3,:] + res = src + res['coord_frame'] = dest + + res['rr'] = np.dot(np.c_[res['rr'], np.ones((res['np'], 1))], t.T) + res['nn'] = np.dot(np.c_[res['nn'], np.zeros((res['np'], 1))], t.T) + return res + + +def transform_meg_chs(chs, trans): + """ + % + % [res, count] = fiff_transform_meg_chs(chs,trans) + % + % Move to another coordinate system in MEG channel channel info + % Count gives the number of channels transformed + % + % NOTE: Only the coil_trans field is modified by this routine, not + % loc which remains to reflect the original data read from the fif file + % + % + + XXX + """ + + res = copy.copy(chs) + + count = 0 + t = trans['trans'] + for ch in res: + if (ch['kind'] == FIFF.FIFFV_MEG_CH + or ch['kind'] == FIFF.FIFFV_REF_MEG_CH): + if (ch['coord_frame'] == trans['from_'] + and ch['coil_trans'] is not None): + ch['coil_trans'] = np.dot(t, ch['coil_trans']) + ch['coord_frame'] = trans['to'] + count += 1 + + if count > 0: + print '\t%d MEG channel locations transformed' % count + + return res, count + + +def transform_eeg_chs(chs, trans): + """ + % + % [res, count] = fiff_transform_eeg_chs(chs,trans) + % + % Move to another coordinate system in EEG channel channel info + % Count gives the number of channels transformed + % + % NOTE: Only the eeg_loc field is modified by this routine, not + % loc which remains to reflect the original data read from the fif file + % + + XXX + """ + res = copy.copy(chs) + + count = 0 + # + # Output unaugmented vectors from the transformation + # + t = trans['trans'][:3,:] + for ch in res: + if ch['kind'] == FIFF.FIFFV_EEG_CH: + if (ch['coord_frame'] == trans['from_'] + and ch['eeg_loc'] is not None): + # + # Transform the augmented EEG location vectors + # + for p in range(ch['eeg_loc'].shape[1]): + ch['eeg_loc'][:, p] = np.dot(t, + np.r_[ch['eeg_loc'][:,p], 1]) + count += 1 + ch['coord_frame'] = trans['to'] + + if count > 0: + print '\t%d EEG electrode locations transformed\n' % count + + return res, count -- 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
