[med-svn] [python-mne] 75/353: apply_forward returns Evoked, cleaned up code
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 6267e6f8f1aedbc3847c38acfffd51d63eb3117e
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date: Mon Feb 6 16:51:17 2012 -0500
apply_forward returns Evoked, cleaned up code
---
mne/forward.py | 134 +++++++++++++++++++++++++++++-----------------
mne/tests/test_forward.py | 50 ++++++++++++++---
2 files changed, 129 insertions(+), 55 deletions(-)
diff --git a/mne/forward.py b/mne/forward.py
index 986c7b7..6094255 100644
--- a/mne/forward.py
+++ b/mne/forward.py
@@ -451,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 source space 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)
@@ -465,16 +465,75 @@ def _stc_src_sel(src, stc):
return src_sel
-def apply_forward(fwd, stc, start=None, stop=None, include=[], exclude=[]):
+def _fill_measurement_info(info, fwd, ch_names, sfreq):
+ """ Fill the measurement info of a Raw or Evoked object
+ """
+ info['bads'] = []
+ info['ch_names'] = ch_names
+ info['chs'] = [deepcopy(ch) for ch in fwd['chs'] if ch['ch_name'] in
+ ch_names]
+ info['nchan'] = len(info['chs'])
+
+ info['filename'] = None
+ info['meas_id'] = None #XXX is this the right thing to do?
+ info['file_id'] = None #XXX is this the right thing to do?
+
+ now = time()
+ sec = np.floor(now)
+ usec = 1e6 * (now - sec)
+
+ info['meas_date'] = np.array([sec, usec], dtype=np.int32)
+ info['highpass'] = np.array(0.0, dtype=np.float32)
+ info['lowpass'] = np.array(sfreq / 2.0, dtype=np.float32)
+ info['sfreq'] = np.array(sfreq, dtype=np.float32)
+ info['projs'] = []
+
+
+def _apply_forward(fwd, stc, start=None, stop=None, include=[], exclude=[]):
+ """ Apply forward model and return data, times, ch_names
+ """
+ if fwd['source_ori'] != FIFF.FIFFV_MNE_FIXED_ORI:
+ raise ValueError('Only fixed-orientation forward operators are '
+ 'supported')
+
+ if np.all(stc.data > 0):
+ warnings.warn('Source estimate only contains currents with positive '
+ 'values. Use pick_normal=True when computing the '
+ 'inverse to compute currents not current magnitudes.')
+
+ fwd = pick_channels_forward(fwd, include=include, exclude=exclude)
+
+ src_sel = _stc_src_sel(fwd['src'], stc)
+
+ gain = fwd['sol']['data'][:, src_sel]
+
+ print 'Projecting source estimate to sensor space...',
+ data = np.dot(gain, stc.data[:, start:stop])
+ print '[done]'
+
+ times = deepcopy(stc.times[start:stop])
+ ch_names = deepcopy(fwd['sol']['row_names'])
+
+ return data, times, ch_names
+
+
+def apply_forward(fwd, stc, evoked_template, start=None, stop=None,
+ include=[], exclude=[]):
"""
Project source space currents to sensor space using a forward operator.
+ The function returns an Evoked object, which is constructed from
+ evoked_template. The evoked_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.
+ evoked_template: Evoked object
+ Evoked 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
@@ -487,40 +546,33 @@ def apply_forward(fwd, stc, start=None, stop=None, include=[], exclude=[]):
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.
+ evoked: Evoked
+ Evoked object with computed sensor space data.
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')
-
- if np.all(stc.data > 0):
- warnings.warn('Source estimate only contains currents with positive '
- 'values. Use pick_normal=True when computing the '
- 'inverse to compute currents not current magnitudes.')
- fwd = pick_channels_forward(fwd, include=include, exclude=exclude)
+ # project the source estimate to the sensor space
+ data, times, ch_names = _apply_forward(fwd, stc, start, stop,
+ include, exclude)
- src_sel = _stc_src_sel(fwd['src'], stc)
+ # store sensor data in an Evoked object using the template
+ evoked = deepcopy(evoked_template)
- gain = fwd['sol']['data'][:, src_sel]
+ evoked.nave = 1
+ evoked.data = data
+ evoked.times = times
- print 'Projecting source estimate to sensor space...',
- data = np.dot(gain, stc.data[:, start:stop])
- print '[done]'
+ sfreq = float(1.0 / stc.tstep)
+ evoked.first = int(np.round(evoked.times[0] * sfreq))
+ evoked.last = evoked.first + evoked.data.shape[1] - 1
- times = deepcopy(stc.times[start:stop])
- ch_names = deepcopy(fwd['sol']['row_names'])
+ # fill the measurement info
+ _fill_measurement_info(evoked.info, fwd, ch_names, sfreq)
- return data, times, ch_names
+ return evoked
def apply_forward_raw(fwd, stc, raw_template, start=None, stop=None,
@@ -557,41 +609,25 @@ def apply_forward_raw(fwd, stc, raw_template, start=None, stop=None,
See Also
--------
- apply_forward: Compute sensor space data and return result as ndarray.
+ apply_forward: Compute sensor space data and return an Evoked object.
"""
# project the source estimate to the sensor space
- data, times, ch_names = apply_forward(fwd, stc, start, stop, include,
- exclude)
+ data, times, ch_names = _apply_forward(fwd, stc, start, stop,
+ include, exclude)
- # store sensor data in Raw object using a template
+ # store sensor data in Raw object using the 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
+ sfreq = float(1.0 / stc.tstep)
raw.first_samp = int(np.round(raw._times[0] * sfreq))
raw.last_samp = raw.first_samp + raw._data.shape[1] - 1
+ # fill the measurement info
+ _fill_measurement_info(raw.info, fwd, ch_names, sfreq)
+
return raw
diff --git a/mne/tests/test_forward.py b/mne/tests/test_forward.py
index 2761adf..0151aac 100644
--- a/mne/tests/test_forward.py
+++ b/mne/tests/test_forward.py
@@ -4,9 +4,10 @@ import numpy as np
from numpy.testing import assert_array_almost_equal, assert_equal
from ..datasets import sample
-from ..fiff import Raw, pick_channels
+from ..fiff import Raw, Evoked, pick_channels
from ..minimum_norm.inverse import _make_stc
-from .. import read_forward_solution, apply_forward_raw, SourceEstimate
+from .. import read_forward_solution, apply_forward, apply_forward_raw,\
+ SourceEstimate
examples_folder = op.join(op.dirname(__file__), '..', '..', 'examples')
@@ -16,6 +17,9 @@ fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg-oct-6-fwd.fif')
fname_raw = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
'test_raw.fif')
+fname_evoked = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data',
+ 'test-ave.fif')
+
def test_io_forward():
"""Test IO for forward solutions
@@ -42,11 +46,45 @@ def test_apply_forward():
stc_data = np.ones((len(vertno[0]) + len(vertno[1]), n_times))
stc = _make_stc(stc_data, t_start, 1.0 / sfreq, vertno)
+ evoked = Evoked(fname_evoked, setno=0)
+
+ evoked = apply_forward(fwd, stc, evoked, start=start, stop=stop)
+
+ data = evoked.data
+ times = evoked.times
+
+ sel = pick_channels(fwd['sol']['row_names'],
+ include=evoked.info['ch_names'])
+
+ gain_sum = np.sum(fwd['sol']['data'][sel, :], axis=1)
+
+ # do some tests
+ assert_array_almost_equal(evoked.info['sfreq'], sfreq)
+ assert_array_almost_equal(np.sum(data, axis=1), n_times * gain_sum)
+ assert_array_almost_equal(times[0], t_start)
+ assert_array_almost_equal(times[-1], t_start + (n_times - 1) / sfreq)
+
+
+def test_apply_forward_raw():
+ """Test projection of source space data to sensor space (Raw)
+ """
+ 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[:, :]
+ data, times = raw_proj[:, :]
sel = pick_channels(fwd['sol']['row_names'],
include=raw_proj.info['ch_names'])
@@ -54,7 +92,7 @@ def test_apply_forward():
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)
+ assert_array_almost_equal(np.sum(data, axis=1), n_times * gain_sum)
+ assert_array_almost_equal(times[0], t_start)
+ assert_array_almost_equal(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