[med-svn] [python-mne] 72/353: ENH: apply_forward, finished implementation, added test
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:33 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit 959221067a2594e2327a30abca384b87743bfb89
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date: Thu Feb 2 14:49:11 2012 -0500
ENH: apply_forward, finished implementation, added test
---
doc/source/whats_new.rst | 2 +
mne/__init__.py | 2 +-
mne/forward.py | 120 ++++++++++++++++++++++++++++++++++++++++++++--
mne/tests/test_forward.py | 43 ++++++++++++++++-
4 files changed, 160 insertions(+), 7 deletions(-)
diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst
index 08814a2..0a67d14 100644
--- a/doc/source/whats_new.rst
+++ b/doc/source/whats_new.rst
@@ -17,6 +17,8 @@ Changelog
- Support of arithmetic of Evoked data (useful to concatenate between runs and compute contrasts) by `Alex Gramfort`_.
+ - Support for computing sensor space data from a source estimate using an MNE forward solution by `Martin Luessi`_.
+
Version 0.2
-----------
diff --git a/mne/__init__.py b/mne/__init__.py
index df79675..6dbb5b1 100644
--- a/mne/__init__.py
+++ b/mne/__init__.py
@@ -4,7 +4,7 @@ from .cov import read_cov, write_cov, write_cov_file, Covariance, \
compute_raw_data_covariance, compute_covariance
from .event import read_events, write_events, find_events, merge_events, \
pick_events
-from .forward import read_forward_solution, apply_forward
+from .forward import read_forward_solution, apply_forward, apply_forward_raw
from .source_estimate import read_stc, write_stc, read_w, write_w, \
SourceEstimate, morph_data, \
spatio_temporal_src_connectivity, \
diff --git a/mne/forward.py b/mne/forward.py
index 85bbd42..986c7b7 100644
--- a/mne/forward.py
+++ b/mne/forward.py
@@ -1,11 +1,15 @@
# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
# Matti Hamalainen <msh at nmr.mgh.harvard.edu>
+# Martin Luessi <mluessi at nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)
+from time import time
+import warnings
+from copy import deepcopy
+
import numpy as np
from scipy import linalg
-import warnings
from .fiff.constants import FIFF
from .fiff.open import fiff_open
@@ -447,7 +451,7 @@ def compute_depth_prior_fixed(G, exp=0.8, limit=10.0):
def _stc_src_sel(src, stc):
- """Select the vertex indices of a forward solution using a source estimate
+ """Select the vertex indices of a source space using a source estimate
"""
src_sel_lh = np.intersect1d(src[0]['vertno'], stc.vertno[0])
src_sel_lh = np.searchsorted(src[0]['vertno'], src_sel_lh)
@@ -462,7 +466,38 @@ def _stc_src_sel(src, stc):
def apply_forward(fwd, stc, start=None, stop=None, include=[], exclude=[]):
+ """
+ Project source space currents to sensor space using a forward operator.
+ Parameters
+ ----------
+ forward: dict
+ Forward operator to use. Has to be fixed-orientation.
+ stc: SourceEstimate
+ The source estimate from which the sensor space data is computed.
+ start: int, optional
+ Index of first time sample (index not time is seconds).
+ stop: int, optional
+ Index of first time sample not to include (index not time is seconds).
+ include: list, optional
+ List of names of channels to include in output. If empty all channels
+ are included.
+ exclude: list, optional
+ List of names of channels to exclude. If empty include all channels.
+
+ Returns
+ -------
+ data: ndarray (n_channels x n_times)
+ Sensor space data.
+ times: ndarray (n_times)
+ Time points in seconds.
+ ch_names: list
+ List with channel names.
+
+ See Also
+ --------
+ apply_forward_raw: Compute sensor space data and return a Raw object.
+ """
if fwd['source_ori'] != FIFF.FIFFV_MNE_FIXED_ORI:
raise ValueError('Only fixed-orientation forward operators are '
'supported')
@@ -479,7 +514,84 @@ def apply_forward(fwd, stc, start=None, stop=None, include=[], exclude=[]):
gain = fwd['sol']['data'][:, src_sel]
print 'Projecting source estimate to sensor space...',
- sens_data = np.dot(gain, stc.data[:, start:stop])
+ data = np.dot(gain, stc.data[:, start:stop])
print '[done]'
- return sens_data
+ times = deepcopy(stc.times[start:stop])
+ ch_names = deepcopy(fwd['sol']['row_names'])
+
+ return data, times, ch_names
+
+
+def apply_forward_raw(fwd, stc, raw_template, start=None, stop=None,
+ include=[], exclude=[]):
+ """
+ Project source space currents to sensor space using a forward operator.
+
+ The function returns a Raw object, which is constructed from raw_template.
+ The raw_template should be from the same MEG system on which the original
+ data was acquired.
+
+ Parameters
+ ----------
+ forward: dict
+ Forward operator to use. Has to be fixed-orientation.
+ stc: SourceEstimate
+ The source estimate from which the sensor space data is computed.
+ raw_template: Raw object
+ Raw object used as template to generate the output argument.
+ start: int, optional
+ Index of first time sample (index not time is seconds).
+ stop: int, optional
+ Index of first time sample not to include (index not time is seconds).
+ include: list, optional
+ List of names of channels to include in output. If empty all channels
+ are included.
+ exclude: list, optional
+ List of names of channels to exclude. If empty include all channels.
+
+ Returns
+ -------
+ raw: Raw object
+ Raw object with computed sensor space data.
+
+ See Also
+ --------
+ apply_forward: Compute sensor space data and return result as ndarray.
+ """
+
+ # project the source estimate to the sensor space
+ data, times, ch_names = apply_forward(fwd, stc, start, stop, include,
+ exclude)
+
+ # store sensor data in Raw object using a template
+ raw = deepcopy(raw_template)
+
+ sfreq = float(1.0 / stc.tstep)
+ now = time()
+ sec = np.floor(now)
+ usec = 1e6 * (now - sec)
+
+ raw.info['bads'] = []
+ raw.info['ch_names'] = ch_names
+ raw.info['chs'] = [deepcopy(ch) for ch in fwd['chs'] if ch['ch_name'] in
+ ch_names]
+ raw.info['nchan'] = len(raw.info['chs'])
+
+ raw.info['filename'] = None
+ raw.info['meas_id'] = None #XXX is this the right thing to do?
+ raw.info['file_id'] = None #XXX is this the right thing to do?
+ raw.info['meas_date'] = np.array([sec, usec], dtype=np.int32)
+ raw.info['highpass'] = np.array(0.0, dtype=np.float32)
+ raw.info['lowpass'] = np.array(sfreq / 2.0, dtype=np.float32)
+ raw.info['sfreq'] = np.array(sfreq, dtype=np.float32)
+ raw.info['projs'] = []
+
+ raw._preloaded = True
+ raw._data = data
+ raw._times = times
+
+ raw.first_samp = int(np.round(raw._times[0] * sfreq))
+ raw.last_samp = raw.first_samp + raw._data.shape[1] - 1
+
+ return raw
diff --git a/mne/tests/test_forward.py b/mne/tests/test_forward.py
index 8fb9222..a825f4b 100644
--- a/mne/tests/test_forward.py
+++ b/mne/tests/test_forward.py
@@ -1,14 +1,20 @@
import os.path as op
-# from numpy.testing import assert_array_almost_equal, assert_equal
+import numpy as np
+from numpy.testing import assert_array_almost_equal, assert_equal
from ..datasets import sample
-from .. import read_forward_solution
+from ..fiff import Raw, pick_channels
+from ..minimum_norm.inverse import _make_stc
+from .. import read_forward_solution, apply_forward_raw, SourceEstimate
+
examples_folder = op.join(op.dirname(__file__), '..', '..', 'examples')
data_path = sample.data_path(examples_folder)
fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg-oct-6-fwd.fif')
+fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
+
def test_io_forward():
"""Test IO for forward solutions
@@ -18,3 +24,36 @@ def test_io_forward():
fwd = read_forward_solution(fname, surf_ori=True)
leadfield = fwd['sol']['data']
# XXX : test something
+
+
+def test_apply_forward():
+ """Test projection of source space data to sensor space
+ """
+ start = 0
+ stop = 5
+ n_times = stop - start - 1
+ sfreq = 10.0
+ t_start = 0.123
+
+ fwd = read_forward_solution(fname, force_fixed=True)
+
+ vertno = [fwd['src'][0]['vertno'], fwd['src'][1]['vertno']]
+ stc_data = np.ones((len(vertno[0]) + len(vertno[1]), n_times))
+ stc = _make_stc(stc_data, t_start, 1.0 / sfreq, vertno)
+
+ raw = Raw(fname_raw)
+
+ raw_proj = apply_forward_raw(fwd, stc, raw, start=start, stop=stop)
+
+ proj_data, proj_times = raw_proj[:, :]
+
+ sel = pick_channels(fwd['sol']['row_names'],
+ include=raw_proj.info['ch_names'])
+
+ gain_sum = np.sum(fwd['sol']['data'][sel, :], axis=1)
+
+ # do some tests
+ assert_array_almost_equal(np.sum(proj_data, axis=1), n_times * gain_sum)
+ assert_array_almost_equal(raw_proj.info['sfreq'], sfreq)
+ assert_array_almost_equal(proj_times[0], t_start)
+ assert_array_almost_equal(proj_times[-1], t_start + (n_times - 1) / sfreq)
--
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