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 a480c5a802f25605b510b2b6fdeb2eaa40c62306 Author: Alexandre Gramfort <[email protected]> Date: Wed Dec 29 10:04:27 2010 -0500 first working example for reading forward solutions --- examples/read_forward.py | 32 +++ fiff/__init__.py | 1 + fiff/forward.py | 538 +++++++++++++++++++++++++++++++++++++++++++++++ fiff/source_space.py | 254 ++++++++++++++++++++++ 4 files changed, 825 insertions(+) diff --git a/examples/read_forward.py b/examples/read_forward.py new file mode 100644 index 0000000..28ad715 --- /dev/null +++ b/examples/read_forward.py @@ -0,0 +1,32 @@ +"""Reading a forward operator a.k.a. lead field matrix +""" +print __doc__ + +import fiff + +fname = 'MNE-sample-data/subjects/sample/bem/sample-5120-bem-sol.fif' +fname = 'sm01a5-ave-oct-6-fwd.fif' + +data = fiff.read_forward_solution(fname) +leadfield = data['sol']['data'] + +print "Leadfield size : %d x %d" % leadfield.shape + +############################################################################### +# Show result + +import pylab as pl +pl.matshow(leadfield[:306,:500]) +pl.xlabel('sources') +pl.ylabel('sensors') +pl.title('Lead field matrix') +pl.show() + +# 3D source space +lh_points = data['src'][0]['rr'] +lh_faces = data['src'][0]['use_tris'] +rh_points = data['src'][1]['rr'] +rh_faces = data['src'][1]['use_tris'] +from enthought.mayavi import mlab +mlab.triangular_mesh(lh_points[:,0], lh_points[:,1], lh_points[:,2], lh_faces) +mlab.triangular_mesh(rh_points[:,0], rh_points[:,1], rh_points[:,2], rh_faces) diff --git a/fiff/__init__.py b/fiff/__init__.py index 552760e..3e77549 100644 --- a/fiff/__init__.py +++ b/fiff/__init__.py @@ -4,4 +4,5 @@ from .evoked import read_evoked from .cov import read_cov from .raw import setup_read_raw, read_raw_segment, read_raw_segment_times from .event import read_events +from .forward import read_forward_solution diff --git a/fiff/forward.py b/fiff/forward.py new file mode 100644 index 0000000..dd6bbdf --- /dev/null +++ b/fiff/forward.py @@ -0,0 +1,538 @@ +import copy +import numpy as np +from scipy import linalg + +from .constants import FIFF +from .open import fiff_open +from .tree import dir_tree_find +from .channels import read_bad_channels +from .tag import find_tag, has_tag +from .source_space import read_source_spaces + + +def transpose_named_matrix(mat): + """Transpose mat inplace (no copy) + """ + mat['nrow'] = mat['ncol'] + mat['ncol'] = mat['nrow'] + mat['row_names'] = mat['col_names'] + mat['col_names'] = mat['row_names'] + mat['data'] = mat['data'].T + return mat + + +def block_diag(A, n): + """ + % + % function bd = mne_block_diag(A,n) + % + % Make or extract a sparse block diagonal matrix + % + % If A is not sparse, then returns a sparse block diagonal "bd", diagonalized from the + % elements in "A". + % "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices. + % Each submatrix is ma x "n", and these submatrices are + % placed down the diagonal of the matrix. + % + % If A is already sparse, then the operation is reversed, yielding a block + % row matrix, where each set of n columns corresponds to a block element + % from the block diagonal. + % + % Routine uses NO for-loops for speed considerations. + """ + from scipy import sparse + + import pdb; pdb.set_trace() + + # if not sparse.issparse(A): # then make block sparse + # ma, na = A.shape + # bdn = na / float(n) # number of submatrices + # + # if bdn > int(bdn): + # raise ValueError, 'Width of matrix must be even multiple of n' + # + # tmp = reshape([1:(ma*bdn)]',ma,bdn); + # i = zeros(ma*n,bdn); + # for iblock = 1:n: + # i((iblock-1)*ma+[1:ma],:) = tmp + # + # i = i.ravel() # row indices foreach sparse bd + # + # j = [1:na]; + # j = j(ones(ma,1),:); + # j = j(:) # column indices foreach sparse bd + # + # bd = sparse(i,j,A(:)); + # else: # already is sparse, unblock it + # [mA,na] = size(A); % matrix always has na columns + # % how many entries in the first column? + # bdn = na/n; % number of blocks + # ma = mA/bdn; % rows in first block + # + # % blocks may themselves contain zero entries. Build indexing as above + # tmp = reshape([1:(ma*bdn)]',ma,bdn); + # i = zeros(ma*n,bdn); + # for iblock = 1:n, + # i((iblock-1)*ma+[1:ma],:) = tmp; + # end + # + # i = i(:); % row indices foreach sparse bd + # + # + # j = [0:mA:(mA*(na-1))]; + # j = j(ones(ma,1),:); + # j = j.ravel() + # + # i += j + # + # bd = full(A(i)); % column vector + # bd = reshape(bd,ma,na); % full matrix + + return bd + + +def transform_source_space_to(src, dest, trans): + """ + % + % [res] = mne_transform_source_space_to(src,dest,trans) + % + % Transform source space data to the desired coordinate system + % + % fname - The name of the file + % include - Include these channels (optional) + % exclude - Exclude these channels (optional) + % + """ + + if src['coord_frame'] == dest: + res = src + return res + + if trans['to'] == src['coord_frame'] and trans['from_'] == dest: + trans = invert_transform(trans) + elif trans['from_'] != src['coord_frame'] or trans['to'] != dest: + raise ValueError, 'Cannot transform the source space using this coordinate transformation' + + t = trans['trans'][:3,:] + res = src + res['coord_frame'] = dest + + res['rr'] = np.dot(np.c_[res['rr'], np.ones((res['np'], 1))], t.T) + res['nn'] = np.dot(np.c_[res['nn'], np.zeros((res['np'], 1))], t.T) + return res + + +def invert_transform(trans): + """ + % + % [itrans] = fiff_invert_transform(trans) + % + % Invert a coordinate transformation + % + """ + itrans = copy.copy(trans) + aux = itrans['from_'] + itrans['from_'] = itrans['to'] + itrans['to'] = aux + itrans['trans'] = linalg.inv(itrans['trans']) + return itrans + + +def read_named_matrix(fid, node, matkind): + """ + % + % [mat] = fiff_read_named_matrix(fid,node) + % + % Read named matrix from the given node + % + """ + + # Descend one level if necessary + if node.block != FIFF.FIFFB_MNE_NAMED_MATRIX: + for k in range(node.nchild): + if node.children[k].block == FIFF.FIFFB_MNE_NAMED_MATRIX: + if has_tag(node.children[k], matkind): + node = node.children[k] + break + else: + raise ValueError, 'Desired named matrix (kind = %d) not available' % matkind + else: + if not has_tag(node, matkind): + raise ValueError, 'Desired named matrix (kind = %d) not available' % matkind + + # Read everything we need + tag = find_tag(fid, node, matkind) + if tag is None: + raise ValueError, 'Matrix data missing' + else: + data = tag.data; + + nrow, ncol = data.shape + tag = find_tag(fid, node, FIFF.FIFF_MNE_NROW) + if tag is not None: + if tag.data != nrow: + raise ValueError, 'Number of rows in matrix data and FIFF_MNE_NROW tag do not match' + + tag = find_tag(fid, node, FIFF.FIFF_MNE_NCOL) + if tag is not None: + if tag.data != ncol: + raise ValueError, 'Number of columns in matrix data and FIFF_MNE_NCOL tag do not match' + + tag = find_tag(fid, node, FIFF.FIFF_MNE_ROW_NAMES) + if tag is not None: + row_names = tag.data.split(':') + else: + row_names = [] + + tag = find_tag(fid, node, FIFF.FIFF_MNE_COL_NAMES) + if tag is not None: + col_names = tag.data.split(':') + else: + col_names = [] + + mat = dict(nrow=nrow, ncol=ncol, row_names=row_names, col_names=col_names, + data=data) + return mat + + +def find_source_space_hemi(src): + """ + % + % function mne_find_source_space_hemi(src) + % + % Return the hemisphere id for a source space + % + % src - The source space to investigate + % hemi - Deduced hemisphere id + % + """ + + xave = src['rr'][:,0].sum(); + + if xave < 0: + hemi = int(FIFF.FIFFV_MNE_SURF_LEFT_HEMI) + else: + hemi = int(FIFF.FIFFV_MNE_SURF_RIGHT_HEMI) + + return hemi + + +def read_one(fid, node): + """ + % + % Read all interesting stuff for one forward solution + % + """ + if node is None: + return None + + tag = find_tag(fid, node, FIFF.FIFF_MNE_SOURCE_ORIENTATION) + if tag is None: + fid.close() + raise ValueError, 'Source orientation tag not found' + + one = dict() + one['source_ori'] = tag.data + tag = find_tag(fid, node, FIFF.FIFF_MNE_COORD_FRAME) + if tag is None: + fid.close() + raise ValueError, 'Coordinate frame tag not found' + + one['coord_frame'] = tag.data + tag = find_tag(fid, node, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS) + if tag is None: + fid.close() + raise ValueError, 'Number of sources not found' + + one['nsource'] = tag.data + tag = find_tag(fid, node, FIFF.FIFF_NCHAN) + if tag is None: + fid.close() + raise ValueError, 'Number of channels not found' + + one['nchan'] = tag.data + try: + one['sol'] = read_named_matrix(fid, node, FIFF.FIFF_MNE_FORWARD_SOLUTION) + one['sol'] = transpose_named_matrix(one['sol']) + except Exception as inst: + fid.close() + raise 'Forward solution data not found (%s)' % inst + + try: + one['sol_grad'] = read_named_matrix(fid,node,FIFF.FIFF_MNE_FORWARD_SOLUTION_GRAD) + one['sol_grad'] = transpose_named_matrix(one['sol_grad']) + except Exception as inst: + one['sol_grad'] = None + + if one['sol']['data'].shape[0] != one['nchan'] or \ + (one['sol']['data'].shape[1] != one['nsource'] and + one['sol']['data'].shape[1] != 3*one['nsource']): + fid.close() + raise ValueError, 'Forward solution matrix has wrong dimensions' + + if one['sol_grad'] is not None: + if one['sol_grad']['data'].shape[0] != one['nchan'] or \ + (one['sol_grad']['data'].shape[1] != 3*one['nsource'] and + one['sol_grad']['data'].shape[1] != 3*3*one['nsource']): + fid.close() + raise ValueError, 'Forward solution gradient matrix has wrong dimensions' + + return one + + +def read_forward_solution(fname, force_fixed=False, surf_ori=False, + include=None, exclude=None): + """ + % + % [fwd] = mne_read_forward_solution(fname,force_fixed,surf_ori,include,exclude) + % + % A forward solution from a fif file + % + % fname - The name of the file + % force_fixed - Force fixed source orientation mode? (optional) + % surf_ori - Use surface based source coordinate system? (optional) + % include - Include these channels (optional) + % exclude - Exclude these channels (optional) + % + """ + + # Open the file, create directory + print 'Reading forward solution from %s...' % fname + fid, tree, _ = fiff_open(fname) + + # Find all forward solutions + fwds = dir_tree_find(tree, FIFF.FIFFB_MNE_FORWARD_SOLUTION) + if len(fwds) == 0: + fid.close() + raise ValueError, 'No forward solutions in %s' % fname + + # Parent MRI data + parent_mri = dir_tree_find(tree,FIFF.FIFFB_MNE_PARENT_MRI_FILE) + if len(fwds) == 0: + fid.close() + raise ValueError, 'No parent MRI information in %s' % fname + parent_mri = parent_mri[0] + + try: + src = read_source_spaces(fid, False, tree) + except Exception as inst: + fid.close() + raise ValueError, 'Could not read the source spaces (%s)' % inst + + for s in src: + s['id'] = find_source_space_hemi(s) + + fwd = None + + # Bad channel list + bads = read_bad_channels(fid, tree) + + print '\t%d bad channels read' % len(bads) + + # Locate and read the forward solutions + megnode = None + eegnode = None + for k in range(len(fwds)): + tag = find_tag(fid, fwds[k], FIFF.FIFF_MNE_INCLUDED_METHODS) + if tag is None: + fid.close() + raise ValueError, 'Methods not listed for one of the forward solutions' + + if tag.data == FIFF.FIFFV_MNE_MEG: + megnode = fwds[k] + elif tag.data == FIFF.FIFFV_MNE_EEG: + eegnode = fwds[k] + + megfwd = read_one(fid, megnode) + if megfwd is not None: + if megfwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI: + ori = 'fixed' + else: + ori = 'free' + + print '\tRead MEG forward solution (%d sources, %d channels, %s orientations)' % ( + megfwd['nsource'], megfwd['nchan'], ori) + + eegfwd = read_one(fid, eegnode) + if eegfwd is not None: + if eegfwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI: + ori = 'fixed' + else: + ori = 'free' + + print '\tRead EEG forward solution (%d sources, %d channels, %s orientations)' % ( + eegfwd['nsource'], eegfwd['nchan'], ori) + + # Merge the MEG and EEG solutions together + if megfwd is not None and eegfwd is not None: + if (megfwd['sol']['data'].shape[1] != eegfwd['sol']['data'].shape[1] or + megfwd['source_ori'] != eegfwd['source_ori'] or + megfwd['nsource'] != eegfwd['nsource'] or + megfwd['coord_frame'] != eegfwd['coord_frame']): + fid.close() + raise ValueError, 'The MEG and EEG forward solutions do not match' + + fwd = megfwd + fwd['sol']['data'] = np.r_[fwd['sol']['data'], eegfwd['sol']['data']] + fwd['sol']['nrow'] = fwd['sol']['nrow'] + eegfwd['sol']['nrow']; + + fwd['sol']['row_names'] = fwd['sol']['row_names'] + eegfwd['sol']['row_names'] + if fwd['sol_grad'] is not None: + fwd['sol_grad']['data'] = np.r_[fwd['sol_grad']['data'], eegfwd['sol_grad']['data']] + fwd['sol_grad']['nrow'] = fwd['sol_grad']['nrow'] + eegfwd['sol_grad']['nrow'] + fwd['sol_grad']['row_names'] = fwd['sol_grad']['row_names'] + eegfwd['sol_grad']['row_names'] + + fwd['nchan'] = fwd['nchan'] + eegfwd['nchan'] + print '\tMEG and EEG forward solutions combined' + elif megfwd is not None: + fwd = megfwd; + else: + fwd = eegfwd; + + del megfwd + del eegfwd + + # Get the MRI <-> head coordinate transformation + tag = find_tag(fid, parent_mri, FIFF.FIFF_COORD_TRANS) + if tag is None: + fid.close() + raise ValueError, 'MRI/head coordinate transformation not found' + else: + mri_head_t = tag.data; + if (mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or + mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): + mri_head_t = invert_transform(mri_head_t) + if (mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI + or mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): + fid.close() + raise ValueError, 'MRI/head coordinate transformation not found' + + fid.close() + + fwd['mri_head_t'] = mri_head_t; + + # Transform the source spaces to the correct coordinate frame + # if necessary + + if (fwd['coord_frame'] != FIFF.FIFFV_COORD_MRI and + fwd['coord_frame'] != FIFF.FIFFV_COORD_HEAD): + raise ValueError, 'Only forward solutions computed in MRI or head coordinates are acceptable' + + nuse = 0 + for s in src: + try: + s = transform_source_space_to(s, fwd['coord_frame'], mri_head_t) + except Exception as inst: + raise ValueError, 'Could not transform source space (%s)' % inst + + nuse += s['nuse'] + + if nuse != fwd['nsource']: + raise ValueError, 'Source spaces do not match the forward solution.' + + print '\tSource spaces transformed to the forward solution coordinate frame' + fwd['src'] = src; + + # Handle the source locations and orientations + if fwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI or force_fixed == True: + nuse = 0 + fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) + fwd['source_nn'] = np.zeros((fwd['nsource'], 3)) + for s in src: + fwd['source_rr'][nuse:nuse+s['nuse'],:] = s['rr'][s['vertno'],:] + fwd['source_nn'][nuse:nuse+s['nuse'],:] = s['nn'][s['vertno'],:] + nuse += s['nuse'] + + # Modify the forward solution for fixed source orientations + if fwd['source_ori'] != FIFF.FIFFV_MNE_FIXED_ORI: + print '\tChanging to fixed-orientation forward solution...' + fix_rot = block_diag(fwd['source_nn'].T, 1) + fwd['sol']['data'] *= fix_rot + fwd['sol']['ncol'] = fwd['nsource'] + fwd['source_ori'] = FIFF.FIFFV_MNE_FIXED_ORI + + if fwd['sol_grad'] is not None: + fwd['sol_grad']['data'] = np.dot(fwd['sol_grad']['data'], np.kron(fix_rot, np.eye(3))) + fwd['sol_grad']['ncol'] = 3*fwd['nsource'] + + print '[done]' + else: + if surf_ori: + # Rotate the local source coordinate systems + print '\tConverting to surface-based source orientations...' + nuse = 0 + pp = 1 + fwd['source_rr'] = np.zeros(fwd['nsource'], 3) + for s in src: + fwd['source_rr'][nuse+1:nuse + s['nuse'],:] = s['rr'][s['vertno'],:] + for p in range(s['nuse']): + # Project out the surface normal and compute SVD + nn = s['nn'][s['vertno'][p],:].T + nn = nn[:,None] + [ U, S, V ] = linalg.svd(np.eye(3,3) - nn*nn.T) + # Make sure that ez is in the direction of nn + if np.sum(nn*U[:,2]) < 0: + U *= -1 + fwd['source_nn'][pp:pp+2,:] = U.T + pp += 3 + nuse += s['nuse'] + + surf_rot = block_diag(fwd['source_nn'].T, 3) + fwd['sol']['data'] = np.dot(fwd['sol']['data'], surf_rot) + if fwd['sol_grad'] is not None: + fwd['sol_grad']['data'] = np.dot(fwd['sol_grad']['data'] * np.kron(surf_rot, np.eye(3))) + + print '[done]' + else: + print '\tCartesian source orientations...' + nuse = 0 + fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) + for s in src: + fwd['source_rr'][nuse:nuse+s['nuse'],:] = s['rr'][s['vertno'],:] + nuse += s['nuse'] + + fwd['source_nn'] = np.kron(np.ones((fwd['nsource'], 1)), np.eye(3)); + print '[done]' + + # Do the channel selection + if include is not None or exclude is not None or bads is not None: + # First do the channels to be included + pick = np.ones(fwd['nchan'], dtype=np.bool) + if include is not None: + for p in range(len(fwd['sol']['row_names'])): + if fwd['sol']['row_names'][p] in include: + pick[p] = True + + # Then exclude what needs to be excluded + if exclude is not None: + for p in range(len(fwd['sol']['row_names'])): + if fwd['sol']['row_names'][p] in exclude: + pick[p] = False + + if bads is not None: + for k in range(len(bads)): + for p in range(len(fwd['sol']['row_names'])): + if fwd['sol']['row_names'][p] in bads: + pick[p] = False + + # Do we have something? + nuse = pick.sum(); + if nuse == 0: + raise ValueError, 'Nothing remains after picking' + + print '\t%d out of %d channels remain after picking\n' % (nuse, fwd['nchan']) + + # Pick the correct rows of the forward operator + fwd['nchan'] = nuse + fwd['sol']['data'] = fwd['sol']['data'][pick == 1,:] + fwd['sol']['nrow'] = nuse + fwd['sol']['row_names'] = [fwd['sol']['row_names'][k] + for k in range(len(pick)) if pick[k]] + + if fwd['sol_grad'] is not None: + fwd['sol_grad']['data'] = fwd['sol_grad']['data'][pick == 1, :] + fwd['sol_grad']['nrow'] = nuse + fwd['sol_grad']['row_names'] = [fwd['sol_grad']['row_names'][k] + for k in range(len(pick)) if pick[k]] + + return fwd \ No newline at end of file diff --git a/fiff/source_space.py b/fiff/source_space.py new file mode 100644 index 0000000..1b5576c --- /dev/null +++ b/fiff/source_space.py @@ -0,0 +1,254 @@ +from math import sqrt +import numpy as np + +from .constants import FIFF +from .tree import dir_tree_find +from .tag import find_tag +from .open import fiff_open + + +def patch_info(nearest): + """ + % + % [pinfo] = mne_patch_info(nearest) + % + % Generate the patch information from the 'nearest' vector in a source space + % + """ + + if nearest is None: + pinfo = None + return pinfo + + indn = np.argsort(nearest) + nearest_sorted = nearest[indn] + + uniq, firsti = np.unique(nearest_sorted, return_index=True) + uniq, lasti = np.unique(nearest_sorted[::-1], return_index=True) + lasti = nearest.size - lasti + + pinfo = list() + for k in range(len(uniq)): + pinfo.append(indn[firsti[k]:lasti[k]]) + + return pinfo + + +def read_source_spaces(source, add_geom=False, tree=None): + """ + % + % [src] = mne_read_source_spaces(source,add_geom,tree) + % + % Reads source spaces from a fif file + % + % source - The name of the file or an open file id + % add_geom - Add geometry information to the source spaces + % tree - Required if source is an open file id, search for the + % source spaces here + % + """ + + # Open the file, create directory + if isinstance(source, str): + fid, tree, _ = fiff_open(source); + open_here = True + else: + fid = source + open_here = False + + # Find all source spaces + spaces = dir_tree_find(tree, FIFF.FIFFB_MNE_SOURCE_SPACE) + if len(spaces) == 0: + if open_here: + fid.close() + raise ValueError, 'No source spaces found' + + src = list() + for s in spaces: + print '\tReading a source space...' + this = _read_one_source_space(fid, s, open_here) + print '[done]' + if add_geom: + complete_source_space_info(this) + + src.append(this) + + print '\t%d source spaces read\n' % len(spaces) + + if open_here: + fid.close() + + return src + + +def _read_one_source_space(fid, this, open_here): + """ + % + % Read all the interesting stuff + % + """ + FIFF_BEM_SURF_NTRI = 3104 + FIFF_BEM_SURF_TRIANGLES = 3106 + + res = dict() + + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_ID) + if tag is None: + res['id'] = int(FIFF.FIFFV_MNE_SURF_UNKNOWN) + else: + res['id'] = int(tag.data) + + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS); + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Number of vertices not found' + + res['np'] = tag.data + + tag = find_tag(fid, this, FIFF_BEM_SURF_NTRI) + if tag is None: + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NTRI) + if tag is None: + res['ntri'] = 0 + else: + res['ntri'] = int(tag.data) + else: + res['ntri'] = tag.data + + tag = find_tag(fid, this, FIFF.FIFF_MNE_COORD_FRAME) + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Coordinate frame information not found' + + res['coord_frame'] = tag.data + + # Vertices, normals, and triangles + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_POINTS) + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Vertex data not found' + + res['rr'] = tag.data + if res['rr'].shape[0] != res['np']: + if open_here: + fid.close() + raise ValueError, 'Vertex information is incorrect' + + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS) + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Vertex normals not found' + + res['nn'] = tag.data + if res['nn'].shape[0] != res['np']: + if open_here: + fid.close() + raise ValueError, 'Vertex normal information is incorrect' + + if res['ntri'] > 0: + tag = find_tag(fid, this, FIFF_BEM_SURF_TRIANGLES) + if tag is None: + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_TRIANGLES) + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Triangulation not found' + else: + res['tris'] = tag.data + else: + res['tris'] = tag.data + + if res['tris'].shape[0] != res['ntri']: + if open_here: + fid.close() + raise ValueError, 'Triangulation information is incorrect' + else: + res['tris'] = None + + # Which vertices are active + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NUSE) + if tag is None: + res['nuse'] = 0; + res['inuse'] = np.zeros(res['nuse'], dtype=np.int) + res['vertno'] = None + else: + res['nuse'] = tag.data + tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_SELECTION) + if tag is None: + if open_here: + fid.close() + raise ValueError, 'Source selection information missing' + + res['inuse'] = tag.data.astype(np.int).T + if len(res['inuse']) != res['np']: + if open_here: + fid.close() + raise ValueError, 'Incorrect number of entries in source space selection' + + res['vertno'] = np.where(res['inuse'])[0] + + # Use triangulation + tag1 = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NUSE_TRI) + tag2 = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES) + if tag1 is None or tag2 is None: + res['nuse_tri'] = 0 + res['use_tris'] = None + else: + res['nuse_tri'] = tag1.data + res['use_tris'] = tag2.data - 1 # index start at 0 in Python + + # Patch-related information + tag1 = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NEAREST) + tag2 = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NEAREST_DIST) + + if tag1 is None or tag2 is None: + res['nearest'] = None + res['nearest_dist'] = None + else: + res['nearest'] = tag1.data + res['nearest_dist'] = tag2.data.T + + res['pinfo'] = patch_info(res['nearest']) + if res['pinfo'] is not None: + print 'Patch information added...' + + return res + + +def complete_source_space_info(this): + """ + """ + # Main triangulation + print '\tCompleting triangulation info...' + this['tri_area'] = np.zeros(this['ntri']) + r1 = this['rr'][this['tris'][:,0],:] + r2 = this['rr'][this['tris'][:,1],:] + r3 = this['rr'][this['tris'][:,2],:] + this['tri_cent'] = (r1 + r2 + r3) / 3.0 + this['tri_nn'] = np.cross((r2-r1), (r3-r1)) + + for p in range(this['ntri']): # XXX : can do better + size = sqrt(np.sum(this['tri_nn'][p,:] * this['tri_nn'][p,:])) + this['tri_area'][p] = size / 2.0 + this['tri_nn'][p,:] = this['tri_nn'][p,:] / size + + print '[done]' + + # Selected triangles + print '\tCompleting selection triangulation info...' + if this['nuse_tri'] > 0: + r1 = this['rr'][this['use_tris'][:,0],:] + r2 = this['rr'][this['use_tris'][:,1],:] + r3 = this['rr'][this['use_tris'][:,2],:] + this['use_tri_cent'] = (r1 + r2 + r3) / 3.0 + this['use_tri_nn'] = np.cross((r2-r1), (r3-r1)); + for p in range(this['nuse_tri']): # XXX can do better + this['use_tri_area'][p] = sqrt(np.sum(this['use_tri_nn'][p,:] * this['use_tri_nn'][p,:])) / 2.0 + + print '[done]' + + -- 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
