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 602fc1d5a9e7dab2ad1c50b98c13dd8a3d75d1dd Author: Alexandre Gramfort <[email protected]> Date: Mon Apr 30 19:53:37 2012 +0200 FIX : fix writing of inv op for mne_analyze + refactoring --- examples/inverse/plot_make_inverse_operator.py | 8 +- mne/fiff/channels.py | 2 +- mne/fiff/cov.py | 4 +- mne/fiff/meas_info.py | 6 +- mne/fiff/pick.py | 15 ++-- mne/forward.py | 119 ++++++++++++++++++++----- mne/minimum_norm/inverse.py | 47 ++++++---- mne/minimum_norm/tests/test_inverse.py | 5 ++ mne/tests/test_forward.py | 4 + 9 files changed, 162 insertions(+), 48 deletions(-) diff --git a/examples/inverse/plot_make_inverse_operator.py b/examples/inverse/plot_make_inverse_operator.py index 6c2a498..c8b55f4 100644 --- a/examples/inverse/plot_make_inverse_operator.py +++ b/examples/inverse/plot_make_inverse_operator.py @@ -19,7 +19,8 @@ import pylab as pl import mne from mne.datasets import sample from mne.fiff import Evoked -from mne.minimum_norm import make_inverse_operator, apply_inverse +from mne.minimum_norm import make_inverse_operator, apply_inverse, \ + write_inverse_operator data_path = sample.data_path('..') fname_fwd = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif' @@ -34,7 +35,7 @@ dSPM = True # Load data evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0)) forward = mne.read_forward_solution(fname_fwd, surf_ori=True) -noise_cov = mne.Covariance(fname_cov) +noise_cov = mne.read_cov(fname_cov) # regularize noise covariance noise_cov = mne.cov.regularize(noise_cov, evoked.info, @@ -43,6 +44,9 @@ noise_cov = mne.cov.regularize(noise_cov, evoked.info, inverse_operator = make_inverse_operator(evoked.info, forward, noise_cov, loose=0.2, depth=0.8) +# Save inverse operator to vizualize with mne_analyze +write_inverse_operator('sample_audvis-eeg-oct-6-eeg-inv.fif', inverse_operator) + # Compute inverse solution stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM, pick_normal=False) diff --git a/mne/fiff/channels.py b/mne/fiff/channels.py index 9e3c708..9ff58c9 100644 --- a/mne/fiff/channels.py +++ b/mne/fiff/channels.py @@ -8,7 +8,7 @@ from .tag import find_tag from .constants import FIFF -def _read_bad_channels(fid, node): +def read_bad_channels(fid, node): """Read bad channels Parameters diff --git a/mne/fiff/cov.py b/mne/fiff/cov.py index 984fcde..3ad01bd 100644 --- a/mne/fiff/cov.py +++ b/mne/fiff/cov.py @@ -11,7 +11,7 @@ from .write import start_block, end_block, write_int, write_name_list, \ from .tag import find_tag from .tree import dir_tree_find from .proj import read_proj, write_proj -from .channels import _read_bad_channels +from .channels import read_bad_channels def read_cov(fid, node, cov_kind): @@ -110,7 +110,7 @@ def read_cov(fid, node, cov_kind): projs = read_proj(fid, this) # Read the bad channel list - bads = _read_bad_channels(fid, this) + bads = read_bad_channels(fid, this) # Put it together cov = dict(kind=cov_kind, diag=diagmat, dim=dim, names=names, diff --git a/mne/fiff/meas_info.py b/mne/fiff/meas_info.py index a85de95..7ebc792 100644 --- a/mne/fiff/meas_info.py +++ b/mne/fiff/meas_info.py @@ -4,6 +4,7 @@ # License: BSD (3-clause) from warnings import warn +from copy import deepcopy import numpy as np from scipy import linalg @@ -13,7 +14,7 @@ from .constants import FIFF from .tag import read_tag from .proj import read_proj, write_proj from .ctf import read_ctf_comp, write_ctf_comp -from .channels import _read_bad_channels +from .channels import read_bad_channels from .write import start_block, end_block, write_string, write_dig_point, \ write_float, write_int, write_coord_trans, write_ch_info, \ @@ -168,7 +169,7 @@ def read_meas_info(fid, tree): comps = read_ctf_comp(fid, meas_info, chs) # Load the bad channel list - bads = _read_bad_channels(fid, meas_info) + bads = read_bad_channels(fid, meas_info) # # Put the data together @@ -326,6 +327,7 @@ def write_meas_info(fid, info, data_type=None): # Channel information for k, c in enumerate(info['chs']): # Scan numbers may have been messed up + c = deepcopy(c) c['scanno'] = k + 1 c['range'] = 1.0 write_ch_info(fid, c) diff --git a/mne/fiff/pick.py b/mne/fiff/pick.py index 78ac7bc..51ec6c0 100644 --- a/mne/fiff/pick.py +++ b/mne/fiff/pick.py @@ -336,11 +336,16 @@ def pick_channels_forward(orig, include=[], exclude=[]): fwd['nchan']) # Pick the correct rows of the forward operator - fwd['nchan'] = nuse fwd['sol']['data'] = fwd['sol']['data'][sel, :] fwd['sol']['nrow'] = nuse - fwd['sol']['row_names'] = [fwd['sol']['row_names'][k] for k in sel] - fwd['chs'] = [fwd['chs'][k] for k in sel] + + ch_names = [fwd['sol']['row_names'][k] for k in sel] + fwd['nchan'] = nuse + fwd['sol']['row_names'] = ch_names + + fwd['info']['chs'] = [fwd['info']['chs'][k] for k in sel] + fwd['info']['nchan'] = nuse + fwd['info']['bads'] = [b for b in fwd['info']['bads'] if b in ch_names] if fwd['sol_grad'] is not None: fwd['sol_grad']['data'] = fwd['sol_grad']['data'][sel, :] @@ -376,11 +381,9 @@ def pick_types_forward(orig, meg=True, eeg=False, include=[], exclude=[]): res : dict Forward solution restricted to selected channel types. """ - info = {'ch_names': orig['sol']['row_names'], 'chs': orig['chs'], - 'nchan': orig['nchan']} + info = orig['info'] 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) diff --git a/mne/forward.py b/mne/forward.py index 0f12412..44b2545 100644 --- a/mne/forward.py +++ b/mne/forward.py @@ -14,9 +14,12 @@ from scipy import linalg from .fiff.constants import FIFF from .fiff.open import fiff_open from .fiff.tree import dir_tree_find +from .fiff.channels import read_bad_channels 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, pick_info, pick_channels +from .fiff.write import write_int, start_block, end_block, \ + write_coord_trans, write_ch_info, write_name_list from .source_space import read_source_spaces_from_tree, find_source_space_hemi from .transforms import transform_source_space_to, invert_transform @@ -161,6 +164,62 @@ def _read_one(fid, node): return one +def read_forward_meas_info(tree, fid): + """Read light measurement info from forward operator + + Parameters + ---------- + tree: tree + FIF tree structure + fid: file id + The file id + + Returns + ------- + info : dict + The measurement info + """ + parent_meg = dir_tree_find(tree, FIFF.FIFFB_MNE_PARENT_MEAS_FILE) + if len(parent_meg) == 0: + fid.close() + raise ValueError('No parent MEG information found in operator') + parent_meg = parent_meg[0] + + # Add channel information + info = dict() + chs = list() + for k in range(parent_meg['nent']): + kind = parent_meg['directory'][k].kind + pos = parent_meg['directory'][k].pos + if kind == FIFF.FIFF_CH_INFO: + tag = read_tag(fid, pos) + chs.append(tag.data) + info['chs'] = chs + + info['ch_names'] = [c['ch_name'] for c in chs] + info['nchan'] = len(chs) + + # Get the MEG device <-> head coordinate transformation + tag = find_tag(fid, parent_meg, FIFF.FIFF_COORD_TRANS) + if tag is None: + fid.close() + raise ValueError('MEG/head coordinate transformation not found') + else: + cand = tag.data + if cand['from'] == FIFF.FIFFV_COORD_DEVICE and \ + cand['to'] == FIFF.FIFFV_COORD_HEAD: + info['dev_head_t'] = cand + elif cand['from'] == FIFF.FIFFV_MNE_COORD_CTF_HEAD and \ + cand['to'] == FIFF.FIFFV_COORD_HEAD: + info['ctf_head_t'] = cand + else: + raise ValueError('MEG device/head coordinate transformation not ' + 'found') + + info['bads'] = read_bad_channels(fid, parent_meg) + return info + + def read_forward_solution(fname, force_fixed=False, surf_ori=False, include=[], exclude=[]): """Read a forward solution a.k.a. lead field @@ -208,21 +267,6 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, raise ValueError('No parent MRI information in %s' % fname) parent_mri = parent_mri[0] - # Parent MEG data - parent_meg = dir_tree_find(tree, FIFF.FIFFB_MNE_PARENT_MEAS_FILE) - if len(parent_meg) == 0: - fid.close() - raise ValueError('No parent MEG information in %s' % fname) - parent_meg = parent_meg[0] - - chs = list() - for k in range(parent_meg['nent']): - kind = parent_meg['directory'][k].kind - pos = parent_meg['directory'][k].pos - if kind == FIFF.FIFF_CH_INFO: - tag = read_tag(fid, pos) - chs.append(tag.data) - try: src = read_source_spaces_from_tree(fid, tree, add_geom=False) except Exception as inst: @@ -317,10 +361,14 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fid.close() raise ValueError('MRI/head coordinate transformation not ' 'found') + fwd['mri_head_t'] = mri_head_t - fid.close() + # + # get parent MEG info + # + fwd['info'] = read_forward_meas_info(tree, fid) - fwd['mri_head_t'] = mri_head_t + fid.close() # Transform the source spaces to the correct coordinate frame # if necessary @@ -416,14 +464,45 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd['surf_ori'] = surf_ori - # Add channel information - fwd['chs'] = chs - fwd = pick_channels_forward(fwd, include=include, exclude=exclude) return fwd +def write_forward_meas_info(fid, info): + """Write measurement info stored in forward solution + + Parameters + ---------- + fid : file id + The file id + info : dict + The measurement info + """ + start_block(fid, FIFF.FIFFB_MNE_PARENT_MEAS_FILE) + # write the MRI <-> head coordinate transformation + if 'dev_head_t' in info: + write_coord_trans(fid, info['dev_head_t']) + if 'ctf_head_t' in info: + write_coord_trans(fid, info['ctf_head_t']) + if 'chs' in info: + # Channel information + write_int(fid, FIFF.FIFF_NCHAN, len(info['chs'])) + for k, c in enumerate(info['chs']): + # Scan numbers may have been messed up + c = deepcopy(c) + c['scanno'] = k + 1 + # c['range'] = 1.0 + write_ch_info(fid, c) + if 'bads' in info: + # Bad channels + if len(info['bads']) > 0: + start_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS) + write_name_list(fid, FIFF.FIFF_MNE_CH_NAME_LIST, info['bads']) + end_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS) + end_block(fid, FIFF.FIFFB_MNE_PARENT_MEAS_FILE) + + def compute_depth_prior(G, exp=0.8, limit=10.0): """Compute weighting for depth prior """ diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index a0a1102..1b86ac9 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -21,9 +21,10 @@ from ..fiff.write import write_int, write_float_matrix, start_file, \ write_coord_trans from ..fiff.cov import read_cov, write_cov -from ..fiff.pick import pick_types +from ..fiff.pick import channel_type from ..cov import prepare_noise_cov -from ..forward import compute_depth_prior, compute_depth_prior_fixed +from ..forward import compute_depth_prior, compute_depth_prior_fixed, \ + read_forward_meas_info, write_forward_meas_info from ..source_space import read_source_spaces_from_tree, \ find_source_space_hemi, _get_vertno, \ write_source_spaces @@ -186,7 +187,6 @@ def read_inverse_operator(fname): # # Read the source spaces # - inv['src'] = read_source_spaces_from_tree(fid, tree, add_geom=False) for s in inv['src']: @@ -213,6 +213,11 @@ def read_inverse_operator(fname): inv['mri_head_t'] = mri_head_t # + # get parent MEG info + # + inv['info'] = read_forward_meas_info(tree, fid) + + # # Transform the source spaces to the correct coordinate frame # if necessary # @@ -331,6 +336,11 @@ def write_inverse_operator(fname, inv): end_block(fid, FIFF.FIFFB_MNE_PARENT_MRI_FILE) # + # Parent MEG measurement info + # + write_forward_meas_info(fid, inv['info']) + + # # Write the source spaces # if 'src' in inv: @@ -927,7 +937,7 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): ---------- info: dict The measurement info to specify the channels to include. - Bad channels in info['bads'] are ignored. + Bad channels in info['bads'] are not used. forward: dict Forward operator noise_cov: Covariance @@ -938,8 +948,6 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): depth: None | float in [0, 1] Depth weighting coefficients. If None, no depth weighting is performed. - # XXX : add support for megreg=0.0, eegreg=0.0 - Returns ------- stc: dict @@ -962,7 +970,7 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): if depth is not None and not (0 < depth <= 1): raise ValueError('depth should be a scalar between 0 and 1') - fwd_ch_names = [c['ch_name'] for c in forward['chs']] + fwd_ch_names = [c['ch_name'] for c in forward['info']['chs']] ch_names = [c['ch_name'] for c in info['chs'] if (c['ch_name'] not in info['bads']) and (c['ch_name'] in fwd_ch_names)] @@ -1008,7 +1016,7 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): source_cov = depth_prior.copy() depth_prior = dict(data=depth_prior, kind=FIFF.FIFFV_MNE_DEPTH_PRIOR_COV, - bads=None, diag=True, names=None, eig=None, + bads=[], diag=True, names=[], eig=None, eigvec=None, dim=depth_prior.size, nfree=1, projs=[]) @@ -1022,7 +1030,7 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): source_cov *= orient_prior orient_prior = dict(data=orient_prior, kind=FIFF.FIFFV_MNE_ORIENT_PRIOR_COV, - bads=None, diag=True, names=None, eig=None, + bads=[], diag=True, names=[], eig=None, eigvec=None, dim=orient_prior.size, nfree=1, projs=[]) else: @@ -1040,8 +1048,8 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): source_cov = dict(data=source_cov, dim=source_cov.size, kind=FIFF.FIFFV_MNE_SOURCE_COV, diag=True, - names=None, projs=[], eig=None, eigvec=None, - nfree=1, bads=None) + names=[], projs=[], eig=None, eigvec=None, + nfree=1, bads=[]) # now np.trace(np.dot(gain, gain.T)) == n_nzero # print np.trace(np.dot(gain, gain.T)), n_nzero @@ -1057,10 +1065,15 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): nave = 1.0 # Handle methods - n_meg = len(pick_types(info, meg=True, eeg=False, exclude=info['bads'])) - n_eeg = len(pick_types(info, meg=False, eeg=True, exclude=info['bads'])) - has_meg = n_meg > 0 - has_eeg = n_eeg > 0 + has_meg = False + has_eeg = False + ch_idx = [k for k, c in enumerate(info['chs']) if c['ch_name'] in ch_names] + for idx in ch_idx: + ch_type = channel_type(info, idx) + if ch_type == 'eeg': + has_eeg = True + if (ch_type == 'mag') or (ch_type == 'grad'): + has_meg = True if has_eeg and has_meg: methods = FIFF.FIFFV_MNE_MEG_EEG elif has_meg: @@ -1079,4 +1092,8 @@ def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): source_nn=forward['source_nn'].copy(), src=deepcopy(forward['src']), fmri_prior=None) + inv_info = deepcopy(forward['info']) + inv_info['bads'] = deepcopy(info['bads']) + inv_op['info'] = inv_info + return inv_op diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 79e6c56..540f4e3 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -112,9 +112,14 @@ def test_apply_inverse_operator(): my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov, loose=0.2, depth=0.8) write_inverse_operator('test-inv.fif', my_inv_op) + read_my_inv_op = read_inverse_operator('test-inv.fif') + _compare(my_inv_op, read_my_inv_op) my_stc = apply_inverse(evoked, my_inv_op, lambda2, dSPM) + assert_true('dev_head_t' in my_inv_op['info']) + assert_true('mri_head_t' in my_inv_op) + assert_equal(stc.times, my_stc.times) assert_array_almost_equal(stc.data, my_stc.data, 2) diff --git a/mne/tests/test_forward.py b/mne/tests/test_forward.py index c235a51..319e300 100644 --- a/mne/tests/test_forward.py +++ b/mne/tests/test_forward.py @@ -1,5 +1,6 @@ import os.path as op +from nose.tools import assert_true import numpy as np from numpy.testing import assert_array_almost_equal, assert_equal @@ -33,6 +34,9 @@ def test_io_forward(): leadfield = fwd['sol']['data'] assert_equal(leadfield.shape, (306, 22494 / 3)) assert_equal(len(fwd['sol']['row_names']), 306) + assert_equal(len(fwd['info']['chs']), 306) + assert_true('dev_head_t' in fwd['info']) + assert_true('mri_head_t' in fwd) def test_apply_forward(): -- 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
