This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit 4100ab9d2de54329eab6bbab182705b5f26e6e79 Author: Alexandre Gramfort <[email protected]> Date: Wed Apr 27 21:24:37 2011 -0400 ENH : adding the possibility to compute MNE on raw data + restriction to 1 label --- examples/plot_minimum_norm_raw_data_in_label.py | 54 ++++++++++ mne/label.py | 7 ++ mne/minimum_norm/__init__.py | 3 +- mne/minimum_norm/inverse.py | 132 +++++++++++++++++++++++- 4 files changed, 192 insertions(+), 4 deletions(-) diff --git a/examples/plot_minimum_norm_raw_data_in_label.py b/examples/plot_minimum_norm_raw_data_in_label.py new file mode 100644 index 0000000..9d5342e --- /dev/null +++ b/examples/plot_minimum_norm_raw_data_in_label.py @@ -0,0 +1,54 @@ +""" +============================================= +Compute MNE-dSPM inverse solution on raw data +============================================= + +Compute dSPM inverse solution on raw dataset restricted +to a brain label and stores the solution in stc files for +visualisation. + +""" + +# Author: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import pylab as pl +import mne +from mne.datasets import sample +from mne.fiff import Raw +from mne.minimum_norm import apply_inverse_raw, read_inverse_operator + + +data_path = sample.data_path('.') +fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' +fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif' +label_name = 'Aud-lh' +fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name + +snr = 3.0 +lambda2 = 1.0 / snr ** 2 +dSPM = True + +# Load data +raw = Raw(fname_raw) +inverse_operator = read_inverse_operator(fname_inv) +label = mne.read_label(fname_label) + +start, stop = raw.time_to_index(0, 15) # read the first 15s of data + +# Compute inverse solution +stc = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM, label, + start, stop) + +# Save result in stc files +stc.save('mne_dSPM_raw_inverse_%s' % label_name) + +############################################################################### +# View activation time-series +pl.plot(1e3 * stc.times, stc.data[::100, :].T) +pl.xlabel('time (ms)') +pl.ylabel('dSPM value') +pl.show() diff --git a/mne/label.py b/mne/label.py index a46f9c5..da5fcad 100755 --- a/mne/label.py +++ b/mne/label.py @@ -33,6 +33,13 @@ def read_label(filename): label['vertices'] = np.array(data[0], dtype=np.int32) label['pos'] = 1e-3 * data[1:4].T label['values'] = data[4] + if filename.endswith('lh.label'): + label['hemi'] = 'lh' + elif filename.endswith('rh.label'): + label['hemi'] = 'rh' + else: + raise ValueError('Cannot find which hemisphere it is. File should end' + ' with lh.label or rh.label') fid.close() return label diff --git a/mne/minimum_norm/__init__.py b/mne/minimum_norm/__init__.py index 1f8924b..273ffa8 100644 --- a/mne/minimum_norm/__init__.py +++ b/mne/minimum_norm/__init__.py @@ -1,2 +1,3 @@ -from .inverse import read_inverse_operator, apply_inverse, minimum_norm +from .inverse import read_inverse_operator, apply_inverse, minimum_norm, \ + apply_inverse_raw from .time_frequency import source_induced_power diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index 0495fb7..ab27cbf 100755 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -14,7 +14,7 @@ from ..fiff.tag import find_tag from ..fiff.matrix import _read_named_matrix, _transpose_named_matrix from ..fiff.proj import read_proj, make_projector from ..fiff.tree import dir_tree_find -from ..fiff.pick import pick_channels_evoked +from ..fiff.pick import pick_channels_evoked, pick_channels from ..cov import read_cov from ..source_space import read_source_spaces_from_tree, find_source_space_hemi @@ -433,8 +433,8 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True): Returns ------- - res: dict - Inverse solution + stc: SourceEstimate + The source estimates """ # @@ -501,6 +501,132 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True): return stc +def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, + label=None, start=None, stop=None, nave=1): + """Apply inverse operator to Raw data + + Computes a L2-norm inverse solution + Actual code using these principles might be different because + the inverse operator is often reused across data sets. + + Parameters + ---------- + raw: Raw object + Evoked data + inverse_operator: dict + Inverse operator read with mne.read_inverse_operator + lambda2: float + The regularization parameter + dSPM: bool + do dSPM ? + label: Label + Restricts the source estimates to a given label + start: int + Index of first time sample (index not time is seconds) + stop: int + Index of last time sample (index not time is seconds) + nave: int + Number of averages used to regularize the solution. + Set to 1 on raw data. + Returns + ------- + stc: SourceEstimate + The source estimates + """ + + # + # Set up the inverse according to the parameters + # + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM) + # + # Pick the correct channels from the data + # + sel = pick_channels(raw.ch_names, include=inv['noise_cov']['names']) + print 'Picked %d channels from the data' % len(sel) + print 'Computing inverse...', + # + # Simple matrix multiplication followed by combination of the + # three current components + # + # This does all the data transformations to compute the weights for the + # eigenleads + # + + src = inv['src'] + lh_vertno = src[0]['vertno'] + rh_vertno = src[1]['vertno'] + + data, times = raw[sel, start:stop] + + trans = inv['reginv'][:, None] * reduce(np.dot, + [inv['eigen_fields']['data'], + inv['whitener'], + inv['proj'], + data]) + + eigen_leads = inv['eigen_leads']['data'] + source_cov = inv['source_cov']['data'][:, None] + noise_norm = inv['noisenorm'][:, None] + + if label is not None: + if label['hemi'] == 'lh': + vertno_sel = np.intersect1d(lh_vertno, label['vertices']) + src_sel = np.searchsorted(lh_vertno, vertno_sel) + lh_vertno = vertno_sel + rh_vertno = np.array([]) + elif label['hemi'] == 'rh': + vertno_sel = np.intersect1d(rh_vertno, label['vertices']) + src_sel = np.searchsorted(rh_vertno, vertno_sel) + len(lh_vertno) + lh_vertno = np.array([]) + rh_vertno = vertno_sel + + noise_norm = noise_norm[src_sel] + + if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: + src_sel = 3 * src_sel + src_sel = np.c_[src_sel, src_sel + 1, src_sel + 2] + src_sel = src_sel.ravel() + + eigen_leads = eigen_leads[src_sel] + source_cov = source_cov[src_sel] + + # + # Transformation into current distributions by weighting the eigenleads + # with the weights computed above + # + if inv['eigen_leads_weighted']: + # + # R^0.5 has been already factored in + # + print '(eigenleads already weighted)...', + sol = np.dot(eigen_leads, trans) + else: + # + # R^0.5 has to factored in + # + print '(eigenleads need to be weighted)...', + sol = np.sqrt(source_cov) * np.dot(eigen_leads, trans) + + if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: + print 'combining the current components...', + sol = combine_xyz(sol) + + if dSPM: + print '(dSPM)...', + sol *= noise_norm + + stc = SourceEstimate(None) + stc.data = sol + stc.tmin = float(times[0]) / raw.info['sfreq'] + stc.tstep = 1.0 / raw.info['sfreq'] + stc.lh_vertno = lh_vertno + stc.rh_vertno = rh_vertno + stc._init_times() + print '[done]' + + return stc + + def _xyz2lf(Lf_xyz, normals): """Reorient leadfield to one component matching the normal to the cortex -- 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
