[med-svn] [python-mne] 55/376: ENH : pick by channel type ENH : data whitening with noise covariance matrix ENH : moving transform functions in transforms.py
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:06 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 4f67f889819ef1b593dc178d389c7d5d3122710b
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
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 <gramfort at nmr.mgh.harvard.edu>
+#
+# 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
More information about the debian-med-commit
mailing list