[med-svn] [python-mne] 166/353: FIX : fix writing of inv op for mne_analyze + refactoring
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:51 UTC 2015
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 <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list