Merge pull request #21 from mluessi/apply_forward

ENH: apply_forward; compute sensor space data from source estimate

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
 -----------
 

mne/__init__.py

@@ -4,7 +4,7 @@
                  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
+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, \

mne/fiff/__init__.py

@@ -11,7 +11,8 @@
 from .raw import Raw, read_raw_segment, read_raw_segment_times, \
                  start_writing_raw, write_raw_buffer, finish_writing_raw
 from .pick import pick_types, pick_channels, pick_types_evoked, \
-                  pick_channels_regexp
+                  pick_channels_regexp, pick_channels_forward, \
+                  pick_types_forward
 
 from .compensator import get_current_comp
 from .proj import compute_spatial_vectors

mne/fiff/pick.py

@@ -315,7 +315,7 @@ def pick_channels_forward(orig, include=[], exclude=[]):
     Returns
     -------
     res : dict
-        Evoked data restricted to selected channels. If include and
+        Forward solution restricted to selected channels. If include and
         exclude are None it returns orig without copy.
     """
 
@@ -351,6 +351,39 @@ def pick_channels_forward(orig, include=[], exclude=[]):
     return fwd
 
 
+def pick_types_forward(orig, meg=True, eeg=False, include=[], exclude=[]):
+    """Pick by channel type and names from a forward operator
+
+    Parameters
+    ----------
+    orig : dict
+        A forward solution
+    meg : bool or string
+        If True include all MEG channels. If False include None
+        If string it can be 'mag' or 'grad' to select only gradiometers
+        or magnetometers. It can also be 'ref_meg' to get CTF
+        reference channels.
+    eeg : bool
+        If True include EEG channels
+    include : list of string
+        List of additional channels to include. If empty do not include any.
+
+    exclude : list of string
+        List of channels to exclude. If empty do not include any.
+
+    Returns
+    -------
+    res : dict
+        Forward solution restricted to selected channel types.
+    """
+    info = {'ch_names': orig['sol']['row_names'], 'chs': orig['chs'],
+            'nchan': orig['nchan']}
+    sel = pick_types(info, meg, eeg, include=include, exclude=exclude)
+    include_ch_names = [info['ch_names'][k] for k in sel]
+
+    return pick_channels_forward(orig, include_ch_names)
+
+
 def channel_indices_by_type(info):
     """Get indices of channels by type
     """

mne/forward.py

@@ -1,8 +1,13 @@
 # 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
 
@@ -11,7 +16,7 @@
 from .fiff.tree import dir_tree_find
 from .fiff.tag import find_tag, read_tag
 from .fiff.matrix import _read_named_matrix, _transpose_named_matrix
-from .fiff.pick import pick_channels_forward
+from .fiff.pick import pick_channels_forward, pick_info, pick_channels
 
 from .source_space import read_source_spaces_from_tree, find_source_space_hemi
 from .transforms import transform_source_space_to, invert_transform
@@ -443,3 +448,192 @@ def compute_depth_prior_fixed(G, exp=0.8, limit=10.0):
     wp = np.minimum(w, wmax)
     depth_prior = (wp / wmax) ** exp
     return depth_prior
+
+
+def _stc_src_sel(src, stc):
+    """ 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)
+
+    src_sel_rh = np.intersect1d(src[1]['vertno'], stc.vertno[1])
+    src_sel_rh = np.searchsorted(src[1]['vertno'], src_sel_rh)\
+                 + len(src[0]['vertno'])
+
+    src_sel = np.r_[src_sel_lh, src_sel_rh]
+
+    return src_sel
+
+
+def _fill_measurement_info(info, fwd, sfreq):
+    """ Fill the measurement info of a Raw or Evoked object
+    """
+    sel = pick_channels(info['ch_names'], fwd['sol']['row_names'])
+    info = pick_info(info, sel)
+    info['bads'] = []
+
+    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'] = []
+
+    return info
+
+
+def _apply_forward(fwd, stc, start=None, stop=None):
+    """ 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.')
+
+    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])
+
+    return data, times
+
+
+def apply_forward(fwd, stc, evoked_template, start=None, stop=None):
+    """
+    Project source space currents to sensor space using a forward operator.
+
+    The sensor space data is computed for all channels present in fwd. Use
+    pick_channels_forward or pick_types_forward to restrict the solution to a
+    subset of channels.
+
+    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. An exception will be raised if the
+    forward operator contains channels that are not present in the template.
+
+
+    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
+        Index of first time sample not to include (index not time is seconds).
+
+    Returns
+    -------
+    evoked: Evoked
+        Evoked object with computed sensor space data.
+
+    See Also
+    --------
+    apply_forward_raw: Compute sensor space data and return a Raw object.
+    """
+
+    # make sure evoked_template contains all channels in fwd
+    for ch_name in fwd['sol']['row_names']:
+        if ch_name not in evoked_template.ch_names:
+            raise ValueError('Channel %s of forward operator not present in '
+                             'evoked_template.' % ch_name)
+
+    # project the source estimate to the sensor space
+    data, times = _apply_forward(fwd, stc, start, stop)
+
+    # store sensor data in an Evoked object using the template
+    evoked = deepcopy(evoked_template)
+
+    evoked.nave = 1
+    evoked.data = data
+    evoked.times = times
+
+    sfreq = float(1.0 / stc.tstep)
+    evoked.first = int(np.round(evoked.times[0] * sfreq))
+    evoked.last = evoked.first + evoked.data.shape[1] - 1
+
+    # fill the measurement info
+    evoked.info = _fill_measurement_info(evoked.info, fwd, sfreq)
+
+    return evoked
+
+
+def apply_forward_raw(fwd, stc, raw_template, start=None, stop=None):
+    """
+    Project source space currents to sensor space using a forward operator.
+
+    The sensor space data is computed for all channels present in fwd. Use
+    pick_channels_forward or pick_types_forward to restrict the solution to a
+    subset of channels.
+
+    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. An exception will be raised if the forward operator
+    contains channels that are not present in the template.
+
+    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).
+
+    Returns
+    -------
+    raw: Raw object
+        Raw object with computed sensor space data.
+
+    See Also
+    --------
+    apply_forward: Compute sensor space data and return an Evoked object.
+    """
+
+    # make sure raw_template contains all channels in fwd
+    for ch_name in fwd['sol']['row_names']:
+        if ch_name not in raw_template.ch_names:
+            raise ValueError('Channel %s of forward operator not present in '
+                             'raw_template.' % ch_name)
+
+    # project the source estimate to the sensor space
+    data, times = _apply_forward(fwd, stc, start, stop)
+
+    # store sensor data in Raw object using the template
+    raw = deepcopy(raw_template)
+
+    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
+    raw.info = _fill_measurement_info(raw.info, fwd, sfreq)
+
+    return raw

mne/tests/test_forward.py

@@ -1,14 +1,26 @@
 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, Evoked, pick_types, pick_types_forward, \
+                   pick_channels_forward
+from ..minimum_norm.inverse import _make_stc
+from .. import read_forward_solution, apply_forward, 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(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
@@ -18,3 +30,66 @@ 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)
+    fwd = pick_types_forward(fwd, meg=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)
+
+    evoked = Evoked(fname_evoked, setno=0)
+
+    evoked = apply_forward(fwd, stc, evoked, start=start, stop=stop)
+
+    data = evoked.data
+    times = evoked.times
+
+    gain_sum = np.sum(fwd['sol']['data'], 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)
+    fwd = pick_types_forward(fwd, meg=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)
+
+    data, times = raw_proj[:, :]
+
+    gain_sum = np.sum(fwd['sol']['data'], axis=1)
+
+    # do some tests
+    assert_array_almost_equal(raw_proj.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)