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 <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
