[med-svn] [python-mne] 225/376: ENH : adding the possibility to compute MNE on raw data + restriction to 1 label
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:44 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 4100ab9d2de54329eab6bbab182705b5f26e6e79
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Apr 27 21:24:37 2011 -0400
ENH : adding the possibility to compute MNE on raw data + restriction to 1 label
---
examples/plot_minimum_norm_raw_data_in_label.py | 54 ++++++++++
mne/label.py | 7 ++
mne/minimum_norm/__init__.py | 3 +-
mne/minimum_norm/inverse.py | 132 +++++++++++++++++++++++-
4 files changed, 192 insertions(+), 4 deletions(-)
diff --git a/examples/plot_minimum_norm_raw_data_in_label.py b/examples/plot_minimum_norm_raw_data_in_label.py
new file mode 100644
index 0000000..9d5342e
--- /dev/null
+++ b/examples/plot_minimum_norm_raw_data_in_label.py
@@ -0,0 +1,54 @@
+"""
+=============================================
+Compute MNE-dSPM inverse solution on raw data
+=============================================
+
+Compute dSPM inverse solution on raw dataset restricted
+to a brain label and stores the solution in stc files for
+visualisation.
+
+"""
+
+# Author: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import pylab as pl
+import mne
+from mne.datasets import sample
+from mne.fiff import Raw
+from mne.minimum_norm import apply_inverse_raw, read_inverse_operator
+
+
+data_path = sample.data_path('.')
+fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
+fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+snr = 3.0
+lambda2 = 1.0 / snr ** 2
+dSPM = True
+
+# Load data
+raw = Raw(fname_raw)
+inverse_operator = read_inverse_operator(fname_inv)
+label = mne.read_label(fname_label)
+
+start, stop = raw.time_to_index(0, 15) # read the first 15s of data
+
+# Compute inverse solution
+stc = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM, label,
+ start, stop)
+
+# Save result in stc files
+stc.save('mne_dSPM_raw_inverse_%s' % label_name)
+
+###############################################################################
+# View activation time-series
+pl.plot(1e3 * stc.times, stc.data[::100, :].T)
+pl.xlabel('time (ms)')
+pl.ylabel('dSPM value')
+pl.show()
diff --git a/mne/label.py b/mne/label.py
index a46f9c5..da5fcad 100755
--- a/mne/label.py
+++ b/mne/label.py
@@ -33,6 +33,13 @@ def read_label(filename):
label['vertices'] = np.array(data[0], dtype=np.int32)
label['pos'] = 1e-3 * data[1:4].T
label['values'] = data[4]
+ if filename.endswith('lh.label'):
+ label['hemi'] = 'lh'
+ elif filename.endswith('rh.label'):
+ label['hemi'] = 'rh'
+ else:
+ raise ValueError('Cannot find which hemisphere it is. File should end'
+ ' with lh.label or rh.label')
fid.close()
return label
diff --git a/mne/minimum_norm/__init__.py b/mne/minimum_norm/__init__.py
index 1f8924b..273ffa8 100644
--- a/mne/minimum_norm/__init__.py
+++ b/mne/minimum_norm/__init__.py
@@ -1,2 +1,3 @@
-from .inverse import read_inverse_operator, apply_inverse, minimum_norm
+from .inverse import read_inverse_operator, apply_inverse, minimum_norm, \
+ apply_inverse_raw
from .time_frequency import source_induced_power
diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py
index 0495fb7..ab27cbf 100755
--- a/mne/minimum_norm/inverse.py
+++ b/mne/minimum_norm/inverse.py
@@ -14,7 +14,7 @@ from ..fiff.tag import find_tag
from ..fiff.matrix import _read_named_matrix, _transpose_named_matrix
from ..fiff.proj import read_proj, make_projector
from ..fiff.tree import dir_tree_find
-from ..fiff.pick import pick_channels_evoked
+from ..fiff.pick import pick_channels_evoked, pick_channels
from ..cov import read_cov
from ..source_space import read_source_spaces_from_tree, find_source_space_hemi
@@ -433,8 +433,8 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True):
Returns
-------
- res: dict
- Inverse solution
+ stc: SourceEstimate
+ The source estimates
"""
#
@@ -501,6 +501,132 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True):
return stc
+def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True,
+ label=None, start=None, stop=None, nave=1):
+ """Apply inverse operator to Raw data
+
+ Computes a L2-norm inverse solution
+ Actual code using these principles might be different because
+ the inverse operator is often reused across data sets.
+
+ Parameters
+ ----------
+ raw: Raw object
+ Evoked data
+ inverse_operator: dict
+ Inverse operator read with mne.read_inverse_operator
+ lambda2: float
+ The regularization parameter
+ dSPM: bool
+ do dSPM ?
+ label: Label
+ Restricts the source estimates to a given label
+ start: int
+ Index of first time sample (index not time is seconds)
+ stop: int
+ Index of last time sample (index not time is seconds)
+ nave: int
+ Number of averages used to regularize the solution.
+ Set to 1 on raw data.
+ Returns
+ -------
+ stc: SourceEstimate
+ The source estimates
+ """
+
+ #
+ # Set up the inverse according to the parameters
+ #
+ inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM)
+ #
+ # Pick the correct channels from the data
+ #
+ sel = pick_channels(raw.ch_names, include=inv['noise_cov']['names'])
+ print 'Picked %d channels from the data' % len(sel)
+ print 'Computing inverse...',
+ #
+ # Simple matrix multiplication followed by combination of the
+ # three current components
+ #
+ # This does all the data transformations to compute the weights for the
+ # eigenleads
+ #
+
+ src = inv['src']
+ lh_vertno = src[0]['vertno']
+ rh_vertno = src[1]['vertno']
+
+ data, times = raw[sel, start:stop]
+
+ trans = inv['reginv'][:, None] * reduce(np.dot,
+ [inv['eigen_fields']['data'],
+ inv['whitener'],
+ inv['proj'],
+ data])
+
+ eigen_leads = inv['eigen_leads']['data']
+ source_cov = inv['source_cov']['data'][:, None]
+ noise_norm = inv['noisenorm'][:, None]
+
+ if label is not None:
+ if label['hemi'] == 'lh':
+ vertno_sel = np.intersect1d(lh_vertno, label['vertices'])
+ src_sel = np.searchsorted(lh_vertno, vertno_sel)
+ lh_vertno = vertno_sel
+ rh_vertno = np.array([])
+ elif label['hemi'] == 'rh':
+ vertno_sel = np.intersect1d(rh_vertno, label['vertices'])
+ src_sel = np.searchsorted(rh_vertno, vertno_sel) + len(lh_vertno)
+ lh_vertno = np.array([])
+ rh_vertno = vertno_sel
+
+ noise_norm = noise_norm[src_sel]
+
+ if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI:
+ src_sel = 3 * src_sel
+ src_sel = np.c_[src_sel, src_sel + 1, src_sel + 2]
+ src_sel = src_sel.ravel()
+
+ eigen_leads = eigen_leads[src_sel]
+ source_cov = source_cov[src_sel]
+
+ #
+ # Transformation into current distributions by weighting the eigenleads
+ # with the weights computed above
+ #
+ if inv['eigen_leads_weighted']:
+ #
+ # R^0.5 has been already factored in
+ #
+ print '(eigenleads already weighted)...',
+ sol = np.dot(eigen_leads, trans)
+ else:
+ #
+ # R^0.5 has to factored in
+ #
+ print '(eigenleads need to be weighted)...',
+ sol = np.sqrt(source_cov) * np.dot(eigen_leads, trans)
+
+ if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI:
+ print 'combining the current components...',
+ sol = combine_xyz(sol)
+
+ if dSPM:
+ print '(dSPM)...',
+ sol *= noise_norm
+
+ stc = SourceEstimate(None)
+ stc.data = sol
+ stc.tmin = float(times[0]) / raw.info['sfreq']
+ stc.tstep = 1.0 / raw.info['sfreq']
+ stc.lh_vertno = lh_vertno
+ stc.rh_vertno = rh_vertno
+ stc._init_times()
+ print '[done]'
+
+ return stc
+
+
def _xyz2lf(Lf_xyz, normals):
"""Reorient leadfield to one component matching the normal to the cortex
--
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