[med-svn] [python-mne] 30/52: ENH : first working support for volume source spaces
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:47 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to annotated tag v0.2
in repository python-mne.
commit c4fd097b774bf73c7803755209e6dc01767f09f2
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Fri Oct 7 16:04:22 2011 -0400
ENH : first working support for volume source spaces
---
mne/minimum_norm/inverse.py | 70 +++++++++++----------------
mne/minimum_norm/tests/test_inverse.py | 19 +++++++-
mne/minimum_norm/time_frequency.py | 10 ++--
mne/parallel.py | 2 +-
mne/source_estimate.py | 86 +++++++++++++++++++++++-----------
mne/source_space.py | 19 ++++++++
mne/transforms.py | 4 +-
7 files changed, 131 insertions(+), 79 deletions(-)
diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py
index 6f31da0..e8c7259 100644
--- a/mne/minimum_norm/inverse.py
+++ b/mne/minimum_norm/inverse.py
@@ -19,7 +19,8 @@ from ..fiff.pick import pick_channels
from ..cov import read_cov, prepare_noise_cov
from ..forward import compute_depth_prior
-from ..source_space import read_source_spaces_from_tree, find_source_space_hemi
+from ..source_space import read_source_spaces_from_tree, \
+ find_source_space_hemi, _get_vertno
from ..transforms import invert_transform, transform_source_space_to
from ..source_estimate import SourceEstimate
@@ -417,14 +418,10 @@ def _assemble_kernel(inv, label, dSPM):
noise_norm = inv['noisenorm'][:, None]
src = inv['src']
- lh_vertno = src[0]['vertno']
- if len(src) > 1:
- rh_vertno = src[1]['vertno']
- else:
- rh_vertno = np.array([])
+ vertno = _get_vertno(src)
if label is not None:
- lh_vertno, rh_vertno, src_sel = _get_label_sel(label, inv)
+ vertno, src_sel = _get_label_sel(label, inv)
noise_norm = noise_norm[src_sel]
@@ -460,42 +457,41 @@ def _assemble_kernel(inv, label, dSPM):
if not dSPM:
noise_norm = None
- return K, noise_norm, lh_vertno, rh_vertno
+ return K, noise_norm, vertno
-def _make_stc(sol, tmin, tstep, lh_vertno, rh_vertno):
+def _make_stc(sol, tmin, tstep, vertno):
stc = SourceEstimate(None)
stc.data = sol
stc.tmin = tmin
stc.tstep = tstep
- stc.lh_vertno = lh_vertno
- stc.rh_vertno = rh_vertno
+ stc.vertno = vertno
stc._init_times()
return stc
def _get_label_sel(label, inv):
src = inv['src']
- lh_vertno = src[0]['vertno']
- if len(src) > 1:
- rh_vertno = src[1]['vertno']
- else:
- rh_vertno = np.array([])
+
+ if src[0]['type'] != 'surf':
+ return Exception('Label are only supported with surface source spaces')
+
+ vertno = [src[0]['vertno'], src[1]['vertno']]
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([])
+ vertno_sel = np.intersect1d(vertno[0], label['vertices'])
+ src_sel = np.searchsorted(vertno[0], vertno_sel)
+ vertno[0] = vertno_sel
+ vertno[1] = 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
+ vertno_sel = np.intersect1d(vertno[1], label['vertices'])
+ src_sel = np.searchsorted(vertno[1], vertno_sel) + len(vertno[0])
+ vertno[0] = np.array([])
+ vertno[1] = vertno_sel
else:
raise Exception("Unknown hemisphere type")
- return lh_vertno, rh_vertno, src_sel
+ return vertno, src_sel
def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True,
@@ -539,7 +535,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True,
print 'Picked %d channels from the data' % len(sel)
print 'Computing inverse...',
- K, noise_norm, _, _ = _assemble_kernel(inv, None, dSPM)
+ K, noise_norm, _ = _assemble_kernel(inv, None, dSPM)
sol = np.dot(K, evoked.data[sel]) # apply imaging kernel
sol = _combine_ori(sol, inv, pick_normal)
@@ -549,14 +545,8 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True,
tstep = 1.0 / evoked.info['sfreq']
tmin = float(evoked.first) / evoked.info['sfreq']
- src = inv['src']
- lh_vertno = src[0]['vertno']
- if len(src) > 1:
- rh_vertno = src[1]['vertno']
- else:
- rh_vertno = np.array([])
-
- stc = _make_stc(sol, tmin, tstep, lh_vertno, rh_vertno)
+ vertno = _get_vertno(inv['src'])
+ stc = _make_stc(sol, tmin, tstep, vertno)
print '[done]'
return stc
@@ -613,16 +603,12 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True,
print 'Picked %d channels from the data' % len(sel)
print 'Computing inverse...',
- src = inv['src']
- lh_vertno = src[0]['vertno']
- rh_vertno = src[1]['vertno']
-
data, times = raw[sel, start:stop]
if time_func is not None:
data = time_func(data)
- K, noise_norm, lh_vertno, rh_vertno = _assemble_kernel(inv, label, dSPM)
+ K, noise_norm, vertno = _assemble_kernel(inv, label, dSPM)
sol = np.dot(K, data)
sol = _combine_ori(sol, inv, pick_normal)
@@ -631,7 +617,7 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True,
tmin = float(times[0])
tstep = 1.0 / raw.info['sfreq']
- stc = _make_stc(sol, tmin, tstep, lh_vertno, rh_vertno)
+ stc = _make_stc(sol, tmin, tstep, vertno)
print '[done]'
return stc
@@ -680,7 +666,7 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True,
print 'Picked %d channels from the data' % len(sel)
print 'Computing inverse...',
- K, noise_norm, lh_vertno, rh_vertno = _assemble_kernel(inv, label, dSPM)
+ K, noise_norm, vertno = _assemble_kernel(inv, label, dSPM)
stcs = list()
tstep = 1.0 / epochs.info['sfreq']
@@ -694,7 +680,7 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True,
if noise_norm is not None:
sol *= noise_norm
- stcs.append(_make_stc(sol, tmin, tstep, lh_vertno, rh_vertno))
+ stcs.append(_make_stc(sol, tmin, tstep, vertno))
print '[done]'
diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py
index 8f29bac..f60dea1 100644
--- a/mne/minimum_norm/tests/test_inverse.py
+++ b/mne/minimum_norm/tests/test_inverse.py
@@ -7,6 +7,7 @@ from ...datasets import sample
from ...label import read_label
from ...event import read_events
from ...epochs import Epochs
+from ...source_estimate import SourceEstimate
from ... import fiff, Covariance, read_forward_solution
from ..inverse import apply_inverse, read_inverse_operator, \
apply_inverse_raw, apply_inverse_epochs, \
@@ -17,6 +18,8 @@ data_path = sample.data_path(examples_folder)
fname_inv = op.join(data_path, 'MEG', 'sample',
# 'sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif')
'sample_audvis-meg-oct-6-meg-inv.fif')
+fname_vol_inv = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis-meg-vol-7-meg-inv.fif')
fname_data = op.join(data_path, 'MEG', 'sample',
'sample_audvis-ave.fif')
fname_cov = op.join(data_path, 'MEG', 'sample',
@@ -63,9 +66,21 @@ def test_inverse_operator():
assert_array_almost_equal(stc.data, my_stc.data, 2)
+def test_inverse_operator_volume():
+ """Test MNE inverse computation on volume source space"""
+ evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0))
+ inverse_operator = read_inverse_operator(fname_vol_inv)
+ stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM)
+ stc.save('tmp-vl.stc')
+ stc2 = SourceEstimate('tmp-vl.stc')
+ assert_true(np.all(stc.data > 0))
+ assert_true(np.all(stc.data < 35))
+ assert_array_almost_equal(stc.data, stc2.data)
+ assert_array_almost_equal(stc.times, stc2.times)
+
+
def test_apply_mne_inverse_raw():
- """Test MNE with precomputed inverse operator on Raw
- """
+ """Test MNE with precomputed inverse operator on Raw"""
start = 3
stop = 10
_, times = raw[0, start:stop]
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index 50af6fa..049d63a 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -65,7 +65,7 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None,
frequencies = np.concatenate([np.arange(band[0], band[1] + df / 2.0, df)
for _, band in bands.iteritems()])
- powers, _, lh_vertno, rh_vertno = _source_induced_power(epochs,
+ powers, _, vertno = _source_induced_power(epochs,
inverse_operator, frequencies,
label=label,
lambda2=lambda2, dSPM=dSPM,
@@ -87,7 +87,7 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None,
verbose=True, copy=False)
tstep = float(decim) / Fs
- stc = _make_stc(power, epochs.times[0], tstep, lh_vertno, rh_vertno)
+ stc = _make_stc(power, epochs.times[0], tstep, vertno)
stcs[name] = stc
print '[done]'
@@ -190,7 +190,7 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
# This does all the data transformations to compute the weights for the
# eigenleads
#
- K, noise_norm, lh_vertno, rh_vertno = _assemble_kernel(inv, label, dSPM)
+ K, noise_norm, vertno = _assemble_kernel(inv, label, dSPM)
if pca:
U, s, Vh = linalg.svd(K, full_matrices=False)
@@ -225,7 +225,7 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
if dSPM:
power *= noise_norm.ravel()[:, None, None] ** 2
- return power, plv, lh_vertno, rh_vertno
+ return power, plv, vertno
def source_induced_power(epochs, inverse_operator, frequencies, label=None,
@@ -280,7 +280,7 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None,
n_jobs: int
Number of jobs to run in parallel
"""
- power, plv, lh_vertno, rh_vertno = _source_induced_power(epochs,
+ power, plv, vertno = _source_induced_power(epochs,
inverse_operator, frequencies,
label, lambda2, dSPM, n_cycles, decim,
use_fft, pick_normal=pick_normal, pca=pca,
diff --git a/mne/parallel.py b/mne/parallel.py
index 12e4164..38a0fc3 100644
--- a/mne/parallel.py
+++ b/mne/parallel.py
@@ -30,7 +30,7 @@ def parallel_func(func, n_jobs, verbose=5):
Number of jobs >= 0
"""
try:
- from scikits.learn.externals.joblib import Parallel, delayed
+ from sklearn.externals.joblib import Parallel, delayed
parallel = Parallel(n_jobs, verbose=verbose)
my_func = delayed(func)
diff --git a/mne/source_estimate.py b/mne/source_estimate.py
index 4ec62e2..6d23912 100644
--- a/mne/source_estimate.py
+++ b/mne/source_estimate.py
@@ -114,23 +114,30 @@ class SourceEstimate(object):
The data in source space
times : array of shape [n_times]
The time vector
- lh_vertno : array of shape [n_dipoles in left hemisphere]
- The indices of the dipoles in the left hemisphere
- rh_vertno : array of shape [n_dipoles in right hemisphere]
- The indices of the dipoles in the right hemisphere
+ vertno : list of array of shape [n_dipoles in each source space]
+ The indices of the dipoles in the different source spaces
"""
def __init__(self, fname):
if fname is not None:
- lh = read_stc(fname + '-lh.stc')
- rh = read_stc(fname + '-rh.stc')
- self.data = np.r_[lh['data'], rh['data']]
- assert lh['tmin'] == rh['tmin']
- assert lh['tstep'] == rh['tstep']
- self.tmin = lh['tmin']
- self.tstep = lh['tstep']
- self.times = self.tmin + self.tstep * np.arange(self.data.shape[1])
- self.lh_vertno = lh['vertices']
- self.rh_vertno = rh['vertices']
+ if fname.endswith('-vl.stc'): # volumne source space
+ vl = read_stc(fname)
+ self.data = vl['data']
+ self.tmin = vl['tmin']
+ self.tstep = vl['tstep']
+ self.times = self.tmin + (self.tstep *
+ np.arange(self.data.shape[1]))
+ self.vertno = [vl['vertices']]
+ else: # surface source spaces
+ lh = read_stc(fname + '-lh.stc')
+ rh = read_stc(fname + '-rh.stc')
+ self.data = np.r_[lh['data'], rh['data']]
+ assert lh['tmin'] == rh['tmin']
+ assert lh['tstep'] == rh['tstep']
+ self.tmin = lh['tmin']
+ self.tstep = lh['tstep']
+ self.times = self.tmin + (self.tstep *
+ np.arange(self.data.shape[1]))
+ self.vertno = [lh['vertices'], rh['vertices']]
def _init_times(self):
"""create self.times"""
@@ -138,18 +145,25 @@ class SourceEstimate(object):
def save(self, fname):
"""save to source estimates to file"""
- lh_data = self.data[:len(self.lh_vertno)]
- rh_data = self.data[-len(self.rh_vertno):]
-
- print 'Writing STC to disk...',
- write_stc(fname + '-lh.stc', tmin=self.tmin, tstep=self.tstep,
- vertices=self.lh_vertno, data=lh_data)
- write_stc(fname + '-rh.stc', tmin=self.tmin, tstep=self.tstep,
- vertices=self.rh_vertno, data=rh_data)
+ if self.is_surface():
+ lh_data = self.data[:len(self.lh_vertno)]
+ rh_data = self.data[-len(self.rh_vertno):]
+
+ print 'Writing STC to disk...',
+ write_stc(fname + '-lh.stc', tmin=self.tmin, tstep=self.tstep,
+ vertices=self.lh_vertno, data=lh_data)
+ write_stc(fname + '-rh.stc', tmin=self.tmin, tstep=self.tstep,
+ vertices=self.rh_vertno, data=rh_data)
+ else:
+ print 'Writing STC to disk...',
+ if not fname.endswith('-vl.stc'):
+ fname += '-vl.stc'
+ write_stc(fname, tmin=self.tmin, tstep=self.tstep,
+ vertices=self.vertno[0], data=self.data)
print '[done]'
def __repr__(self):
- s = "%d vertices" % (len(self.lh_vertno) + len(self.rh_vertno))
+ s = "%d vertices" % sum([len(v) for v in self.vertno])
s += ", tmin : %s (ms)" % (1e3 * self.tmin)
s += ", tmax : %s (ms)" % (1e3 * self.times[-1])
s += ", tstep : %s (ms)" % (1e3 * self.tstep)
@@ -166,6 +180,22 @@ class SourceEstimate(object):
self.data = self.data[:, mask]
self.tmin = self.times[0]
+ @property
+ def lh_vertno(self):
+ return self.vertno[0]
+
+ @property
+ def rh_vertno(self):
+ return self.vertno[1]
+
+ def is_surface(self):
+ """Returns True if source estimate is defined over surfaces
+ """
+ if len(self.vertno) == 1:
+ return False
+ else:
+ return True
+
def __add__(self, a):
stc = copy.deepcopy(self)
stc += a
@@ -383,6 +413,10 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
"""
from scipy import sparse
+ if not stc_from.is_surface():
+ raise ValueError('Morphing is only possible with surface source '
+ 'estimates')
+
if subjects_dir is None:
if 'SUBJECTS_DIR' in os.environ:
subjects_dir = os.environ['SUBJECTS_DIR']
@@ -393,7 +427,6 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
surf_path_from = os.path.join(subjects_dir, subject_from, 'surf')
tris.append(read_surface(os.path.join(surf_path_from, 'lh.sphere.reg'))[1])
tris.append(read_surface(os.path.join(surf_path_from, 'rh.sphere.reg'))[1])
- vertices = [stc_from.lh_vertno, stc_from.rh_vertno]
sphere = os.path.join(subjects_dir, subject_to, 'surf', 'lh.sphere.reg')
lhs = read_surface(sphere)[0]
@@ -412,7 +445,7 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
e.data[e.data == 2] = 1
n_vertices = e.shape[0]
e = e + sparse.eye(n_vertices, n_vertices)
- idx_use = vertices[hemi]
+ idx_use = stc_from.vertno[hemi]
n_iter = 100 # max nb of smoothing iterations
for k in range(n_iter):
e_use = e[:, idx_use]
@@ -456,8 +489,7 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
nearest[1, k:k + dr] = np.argmax(dots, axis=1)
stc_to = copy.deepcopy(stc_from)
- stc_to.lh_vertno = nearest[0]
- stc_to.rh_vertno = nearest[1]
+ stc_to.vertno = [nearest[0], nearest[1]]
stc_to.data = np.r_[dmap[0][nearest[0], :], dmap[1][nearest[1], :]]
print '[done]'
diff --git a/mne/source_space.py b/mne/source_space.py
index 5beaa28..87eb3db 100644
--- a/mne/source_space.py
+++ b/mne/source_space.py
@@ -120,6 +120,18 @@ def _read_one_source_space(fid, this):
else:
res['id'] = int(tag.data)
+ tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_TYPE)
+ if tag is None:
+ raise ValueError('Unknown source space type')
+ else:
+ src_type = int(tag.data)
+ if src_type == 1:
+ res['type'] = 'surf'
+ elif src_type == 2:
+ res['type'] = 'vol'
+ else:
+ raise ValueError('Unknown source space type (%d)' % src_type)
+
tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS)
if tag is None:
raise ValueError('Number of vertices not found')
@@ -273,3 +285,10 @@ def find_source_space_hemi(src):
hemi = int(FIFF.FIFFV_MNE_SURF_RIGHT_HEMI)
return hemi
+
+
+def _get_vertno(src):
+ vertno = list()
+ for s in src:
+ vertno.append(s['vertno'])
+ return vertno
diff --git a/mne/transforms.py b/mne/transforms.py
index 8908e5a..a4b169e 100644
--- a/mne/transforms.py
+++ b/mne/transforms.py
@@ -81,8 +81,8 @@ def transform_coordinates(filename, pos, orig, dest):
Example
-------
- >>> transform_coordinates('all-trans.fif', np.eye(3), 'meg', 'fs_tal')
- >>> transform_coordinates('all-trans.fif', np.eye(3), 'mri', 'mni_tal')
+ transform_coordinates('all-trans.fif', np.eye(3), 'meg', 'fs_tal')
+ transform_coordinates('all-trans.fif', np.eye(3), 'mri', 'mni_tal')
"""
# Read the fif file containing all necessary transformations
fid, tree, directory = fiff_open(filename)
--
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