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