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 <[email protected]> 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 <[email protected]> # Matti Hamalainen <[email protected]> +# Martin Luessi <[email protected]> # # 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
