[med-svn] [python-mne] 18/376: now writing evoked data
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:21:57 UTC 2015
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 69291b016bc4f4a6360f3d6e409e978493487d72
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Thu Dec 30 14:00:21 2010 -0500
now writing evoked data
---
examples/read_evoked.py | 6 ++
fiff/__init__.py | 2 +-
fiff/ctf.py | 102 +++++++++++++++------
fiff/evoked.py | 176 +++++++++++++++++++++++++++++++-----
fiff/open.py | 7 --
fiff/tag.py | 1 -
fiff/tests/test_evoked.py | 32 +++++++
fiff/tree.py | 75 +++++++++++++---
fiff/write.py | 222 +++++++++++++++++++++++++++++++++++++++++++---
9 files changed, 545 insertions(+), 78 deletions(-)
diff --git a/examples/read_evoked.py b/examples/read_evoked.py
index 8fe664f..08a3946 100644
--- a/examples/read_evoked.py
+++ b/examples/read_evoked.py
@@ -8,6 +8,12 @@ fname = 'MNE-sample-data/MEG/sample/sample_audvis-ave.fif'
data = fiff.read_evoked(fname)
+fiff.write_evoked('evoked.fif', data)
+data2 = fiff.read_evoked('evoked.fif')
+
+from scipy import linalg
+print linalg.norm(data['evoked']['epochs'] - data2['evoked']['epochs'])
+
###############################################################################
# Show result
diff --git a/fiff/__init__.py b/fiff/__init__.py
index 75a998c..10be9cc 100644
--- a/fiff/__init__.py
+++ b/fiff/__init__.py
@@ -2,7 +2,7 @@ __version__ = '0.1.git'
from .constants import FIFF
from .open import fiff_open
-from .evoked import read_evoked
+from .evoked import read_evoked, write_evoked
from .cov import read_cov, write_cov, write_cov_file
from .raw import setup_read_raw, read_raw_segment, read_raw_segment_times
from .event import read_events
diff --git a/fiff/ctf.py b/fiff/ctf.py
index 2bfb52f..f5c7461 100644
--- a/fiff/ctf.py
+++ b/fiff/ctf.py
@@ -23,42 +23,45 @@ def read_named_matrix(fid, node, matkind):
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;
+ node = node.children(k)
+ break
else:
- raise ValueError, 'Desired named matrix (kind = %d) not available' % matkind
+ raise ValueError, 'Desired named matrix (kind = %d) not' \
+ ' available' % matkind
else:
- if not has_tag(node,matkind):
- raise 'Desired named matrix (kind = %d) not available' % matkind
+ if not has_tag(node, matkind):
+ raise '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'
+ raise ValueError, 'Matrix data missing'
else:
- data = tag.data;
+ 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'
+ 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'
+ 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;
+ row_names = tag.data
else:
row_names = None
tag = find_tag(fid, node, FIFF.FIFF_MNE_COL_NAMES)
if tag is not None:
- col_names = tag.data;
+ col_names = tag.data
else:
col_names = None
@@ -74,7 +77,7 @@ def read_named_matrix(fid, node, matkind):
else:
mat['col_names'] = None
- mat['data'] = data;
+ mat['data'] = data
return mat
@@ -93,12 +96,12 @@ def read_ctf_comp(fid, node, chs):
for node in comps:
# Read the data we need
- mat = read_named_matrix(fid, node, FIFF.FIFF_MNE_CTF_COMP_DATA)
+ mat = read_named_matrix(fid, node, FIFF.FIFF_MNE_CTF_COMP_DATA)
for p in range(node.nent):
kind = node.dir[p].kind
- pos = node.dir[p].pos
+ pos = node.dir[p].pos
if kind == FIFF.FIFF_MNE_CTF_COMP_KIND:
- tag = read_tag(fid,pos)
+ tag = read_tag(fid, pos)
break
else:
raise ValueError, 'Compensation type not found'
@@ -118,15 +121,15 @@ def read_ctf_comp(fid, node, chs):
for p in range(node.nent):
kind = node.dir[p].kind
- pos = node.dir[p].pos
+ pos = node.dir[p].pos
if kind == FIFF.FIFF_MNE_CTF_COMP_CALIBRATED:
- tag = read_tag(fid,pos)
+ tag = read_tag(fid, pos)
calibrated = tag.data
break
else:
calibrated = False
- one['save_calibrated'] = calibrated;
+ one['save_calibrated'] = calibrated
one['rowcals'] = np.ones(1, mat.shape[0])
one['colcals'] = np.ones(1, mat.shape[1])
if not calibrated:
@@ -143,9 +146,11 @@ def read_ctf_comp(fid, node, chs):
for col in range(mat.data.shape[1]):
p = ch_names.count(mat.col_names[col])
if p == 0:
- raise ValueError, 'Channel %s is not available in data' % mat.col_names[col]
+ raise ValueError, 'Channel %s is not available in data' \
+ % mat.col_names[col]
elif p > 1:
- raise ValueError, 'Ambiguous channel %s' % mat.col_names[col]
+ raise ValueError, 'Ambiguous channel %s' % \
+ mat.col_names[col]
col_cals[col] = 1.0 / (chs[p].range * chs[p].cal)
@@ -154,13 +159,16 @@ def read_ctf_comp(fid, node, chs):
for row in range(mat.data.shape[0]):
p = ch_names.count(mat.row_names[row])
if p == 0:
- raise ValueError, 'Channel %s is not available in data', mat.row_names[row]
+ raise ValueError, 'Channel %s is not available in data' \
+ % mat.row_names[row]
elif p > 1:
- raise ValueError, 'Ambiguous channel %s' % mat.row_names[row]
+ raise ValueError, 'Ambiguous channel %s' % \
+ mat.row_names[row]
row_cals[row] = chs[p].range * chs[p].cal
- mat.data = np.dot(np.diag(row_cals), np.dot(mat.data, np.diag(col_cals)))
+ mat.data = np.dot(np.diag(row_cals), np.dot(mat.data,
+ np.diag(col_cals)))
one.rowcals = row_cals
one.colcals = col_cals
@@ -173,3 +181,47 @@ def read_ctf_comp(fid, node, chs):
print '\tRead %d compensation matrices\n' % len(compdata)
return compdata
+
+
+###############################################################################
+# Writing
+
+from scipy import linalg
+from .write import start_block, end_block, write_int, write_named_matrix
+
+
+def write_ctf_comp(fid, comps):
+ """
+ %
+ % fiff_write_ctf_comp(fid,comps)
+ %
+ % Writes the CTF compensation data into a fif file
+ %
+ % fid An open fif file descriptor
+ % comps The compensation data to write
+ %
+ """
+
+ if len(comps) <= 0:
+ return
+
+ # This is very simple in fact
+ start_block(fid, FIFF.FIFFB_MNE_CTF_COMP)
+ for comp in comps:
+ start_block(fid, FIFF.FIFFB_MNE_CTF_COMP_DATA)
+ # Write the compensation kind
+ write_int(fid, FIFF.FIFF_MNE_CTF_COMP_KIND, comp['ctfkind'])
+ write_int(fid, FIFF.FIFF_MNE_CTF_COMP_CALIBRATED,
+ comp['save_calibrated'])
+
+ # Write an uncalibrated or calibrated matrix
+ import pdb; pdb.set_trace()
+
+ comp['data']['data'] = linalg.inv(
+ np.dot(np.diag(comp['rowcals'].ravel())),
+ np.dot(comp.data.data,
+ linalg.inv(np.diag(comp.colcals.ravel()))))
+ write_named_matrix(fid, FIFF.FIFF_MNE_CTF_COMP_DATA, comp['data'])
+ end_block(fid, FIFF.FIFFB_MNE_CTF_COMP_DATA)
+
+ end_block(fid, FIFF.FIFFB_MNE_CTF_COMP)
diff --git a/fiff/evoked.py b/fiff/evoked.py
index 09cd9f3..013c481 100644
--- a/fiff/evoked.py
+++ b/fiff/evoked.py
@@ -16,22 +16,22 @@ def read_evoked(fname, setno=0):
"""
if setno < 0:
- raise ValueError, 'Data set selector must be positive'
+ raise ValueError, 'Data set selector must be positive'
print 'Reading %s ...\n' % fname
- fid, tree, _ = fiff_open(fname);
+ fid, tree, _ = fiff_open(fname)
# Read the measurement info
info, meas = read_meas_info(fid, tree)
info['filename'] = fname
# Locate the data of interest
- processed = dir_tree_find(meas,FIFF.FIFFB_PROCESSED_DATA);
+ processed = dir_tree_find(meas, FIFF.FIFFB_PROCESSED_DATA)
if len(processed) == 0:
fid.close()
raise ValueError, 'Could not find processed data'
- evoked = dir_tree_find(meas,FIFF.FIFFB_EVOKED)
+ evoked = dir_tree_find(meas, FIFF.FIFFB_EVOKED)
if len(evoked) == 0:
fid.close(fid)
raise ValueError, 'Could not find evoked data'
@@ -54,13 +54,13 @@ def read_evoked(fname, setno=0):
is_smsh = np.c_[is_smsh, np.zeros((1, sets[k]['naspect']))]
naspect += sets[k]['naspect']
- saspects = dir_tree_find(evoked[k], FIFF.FIFFB_SMSH_ASPECT)
+ saspects = dir_tree_find(evoked[k], FIFF.FIFFB_SMSH_ASPECT)
nsaspects = len(saspects)
if nsaspects > 0:
sets[k]['naspect'] += nsaspects
sets[k]['naspect'] = [sets[k]['naspect'], saspects] # XXX : potential bug
is_smsh = np.c_[is_smsh, np.ones(1, sets[k]['naspect'])]
- naspect += nsaspects;
+ naspect += nsaspects
print '\t%d evoked data sets containing a total of %d data aspects' \
' in %s\n' % (len(evoked), naspect, fname)
@@ -96,7 +96,7 @@ def read_evoked(fname, setno=0):
comment = None
for k in range(my_evoked.nent):
kind = my_evoked.directory[k].kind
- pos = my_evoked.directory[k].pos
+ pos = my_evoked.directory[k].pos
if kind == FIFF.FIFF_COMMENT:
tag = read_tag(fid, pos)
comment = tag.data
@@ -118,37 +118,39 @@ def read_evoked(fname, setno=0):
q += 1
if comment is None:
- comment = 'No comment'
+ comment = 'No comment'
# Local channel information?
if nchan > 0:
if chs is None:
fid.close()
- raise ValueError, 'Local channel information was not found when it was expected.'
+ raise ValueError, 'Local channel information was not found ' \
+ 'when it was expected.'
if len(chs) != nchan:
fid.close()
- raise ValueError, 'Number of channels and number of channel definitions are different'
+ raise ValueError, 'Number of channels and number of channel ' \
+ 'definitions are different'
info.chs = chs
info.nchan = nchan
- print '\tFound channel information in evoked data. nchan = %d\n' % nchan
+ print '\tFound channel information in evoked data. nchan = %d' % nchan
if sfreq > 0:
info['sfreq'] = sfreq
nsamp = last - first + 1
- print '\tFound the data of interest:\n'
- print '\t\tt = %10.2f ... %10.2f ms (%s)\n' % (
+ print '\tFound the data of interest:'
+ print '\t\tt = %10.2f ... %10.2f ms (%s)' % (
1000*first/info['sfreq'], 1000*last/info['sfreq'], comment)
if info['comps'] is not None:
- print '\t\t%d CTF compensation matrices available\n' % len(info['comps'])
+ print '\t\t%d CTF compensation matrices available' % len(info['comps'])
# Read the data in the aspect block
nave = 1
epoch = []
for k in range(my_aspect.nent):
kind = my_aspect.directory[k].kind
- pos = my_aspect.directory[k].pos
+ pos = my_aspect.directory[k].pos
if kind == FIFF.FIFF_COMMENT:
tag = read_tag(fid, pos)
comment = tag.data
@@ -171,11 +173,11 @@ def read_evoked(fname, setno=0):
'(nepoch = %d nchan = %d)' % (nepoch, info.nchan)
if nepoch == 1:
- # Only one epoch
- all_data = epoch[0].data
- # May need a transpose if the number of channels is one
- if all_data.shape[1] == 1 and info.nchan == 1:
- all_data = all_data.T
+ # Only one epoch
+ all_data = epoch[0].data
+ # May need a transpose if the number of channels is one
+ if all_data.shape[1] == 1 and info.nchan == 1:
+ all_data = all_data.T
else:
# Put the old style epochs together
@@ -204,3 +206,137 @@ def read_evoked(fname, setno=0):
fid.close()
return data
+
+###############################################################################
+# Writing
+
+from .tree import copy_tree
+from .write import start_file, start_block, end_file, end_block, write_id, \
+ write_float, write_int, write_coord_trans, write_ch_info, \
+ write_dig_point, write_name_list, write_string, \
+ write_float_matrix
+from .proj import write_proj
+from .ctf import write_ctf_comp
+
+
+def write_evoked(name, data):
+ """
+ %
+ % function fiff_write_evoked(name,data)
+ %
+ % name filename
+ % data the data structure returned from fiff_read_evoked
+ %
+ %
+ """
+
+ # Create the file and save the essentials
+ fid = start_file(name)
+ start_block(fid, FIFF.FIFFB_MEAS)
+ write_id(fid, FIFF.FIFF_BLOCK_ID)
+ if data['info']['meas_id'] is not None:
+ write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, data['info']['meas_id'])
+
+ # Measurement info
+ start_block(fid, FIFF.FIFFB_MEAS_INFO)
+
+ # Blocks from the original
+ blocks = [FIFF.FIFFB_SUBJECT, FIFF.FIFFB_HPI_MEAS, FIFF.FIFFB_HPI_RESULT,
+ FIFF.FIFFB_ISOTRAK, FIFF.FIFFB_PROCESSING_HISTORY]
+
+ have_hpi_result = False
+ have_isotrak = False
+ if len(blocks) > 0 and 'filename' in data['info'] and \
+ data['info']['filename'] is not None:
+ fid2, tree, _ = fiff_open(data['info']['filename'])
+ for block in blocks:
+ nodes = dir_tree_find(tree, block)
+ copy_tree(fid2, tree['id'], nodes, fid)
+ if block == FIFF.FIFFB_HPI_RESULT and len(nodes) > 0:
+ have_hpi_result = True
+ if block == FIFF.FIFFB_ISOTRAK and len(nodes) > 0:
+ have_isotrak = True
+ fid2.close()
+
+
+ # General
+ write_float(fid, FIFF.FIFF_SFREQ, data['info']['sfreq'])
+ write_float(fid, FIFF.FIFF_HIGHPASS, data['info']['highpass'])
+ write_float(fid, FIFF.FIFF_LOWPASS, data['info']['lowpass'])
+ write_int(fid, FIFF.FIFF_NCHAN, data['info']['nchan'])
+ if data['info']['meas_date'] is not None:
+ write_int(fid, FIFF.FIFF_MEAS_DATE, data['info']['meas_date'])
+
+ # Coordinate transformations if the HPI result block was not there
+ if not have_hpi_result:
+ if data['info']['dev_head_t'] is not None:
+ write_coord_trans(fid, data['info']['dev_head_t'])
+
+ if data['info']['ctf_head_t'] is not None:
+ write_coord_trans(fid, data['info']['ctf_head_t'])
+
+ # Channel information
+ for k in range(data['info']['nchan']):
+ # Scan numbers may have been messed up
+ data['info']['chs'][k]['scanno'] = k
+ write_ch_info(fid, data['info']['chs'][k])
+
+ # Polhemus data
+ if data['info']['dig'] is not None and not have_isotrak:
+ start_block(fid, FIFF.FIFFB_ISOTRAK)
+ for d in data['info']['dig']:
+ write_dig_point(fid, d)
+
+ end_block(fid, FIFF.FIFFB_ISOTRAK)
+
+ # Projectors
+ write_proj(fid, data['info']['projs'])
+
+ # CTF compensation info
+ write_ctf_comp(fid, data['info']['comps'])
+
+ # Bad channels
+ if len(data['info']['bads']) > 0:
+ start_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS)
+ write_name_list(fid, FIFF.FIFF_MNE_CH_NAME_LIST, data['info']['bads'])
+ end_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS)
+
+ end_block(fid, FIFF.FIFFB_MEAS_INFO)
+
+ # One or more evoked data sets
+ start_block(fid, FIFF.FIFFB_PROCESSED_DATA)
+ data_evoked = data['evoked']
+ if not isinstance(data_evoked, list):
+ data_evoked = [data_evoked]
+
+ for evoked in data_evoked:
+ start_block(fid, FIFF.FIFFB_EVOKED)
+
+ # Comment is optional
+ if len(evoked['comment']) > 0:
+ write_string(fid, FIFF.FIFF_COMMENT, evoked['comment'])
+
+ # First and last sample
+ write_int(fid, FIFF.FIFF_FIRST_SAMPLE, evoked['first'])
+ write_int(fid, FIFF.FIFF_LAST_SAMPLE, evoked['last'])
+
+ # The epoch itself
+ start_block(fid, FIFF.FIFFB_ASPECT)
+
+ write_int(fid, FIFF.FIFF_ASPECT_KIND, evoked['aspect_kind'])
+ write_int(fid, FIFF.FIFF_NAVE, evoked['nave'])
+
+ decal = np.zeros((data['info']['nchan'], data['info']['nchan']))
+ for k in range(data['info']['nchan']): # XXX : can be improved
+ decal[k, k] = 1.0 / data['info']['chs'][k]['cal']
+
+ write_float_matrix(fid, FIFF.FIFF_EPOCH,
+ np.dot(decal, evoked['epochs']))
+ end_block(fid, FIFF.FIFFB_ASPECT)
+ end_block(fid, FIFF.FIFFB_EVOKED)
+
+ end_block(fid, FIFF.FIFFB_PROCESSED_DATA)
+
+ end_block(fid, FIFF.FIFFB_MEAS)
+
+ end_file(fid)
diff --git a/fiff/open.py b/fiff/open.py
index 25ca967..5feaa90 100644
--- a/fiff/open.py
+++ b/fiff/open.py
@@ -9,9 +9,7 @@ def fiff_open(fname, verbose=False):
tag = read_tag_info(fid)
- #
# Check that this looks like a fif file
- #
if tag.kind != FIFF.FIFF_FILE_ID:
raise ValueError, 'file does not start with a file id tag'
@@ -26,9 +24,7 @@ def fiff_open(fname, verbose=False):
if tag.kind != FIFF.FIFF_DIR_POINTER:
raise ValueError, 'file does have a directory pointer'
- #
# Read or create the directory tree
- #
if verbose:
print '\tCreating tag directory for %s...' % fname
@@ -50,10 +46,7 @@ def fiff_open(fname, verbose=False):
if verbose:
print '[done]\n'
- #
# Back to the beginning
- #
fid.seek(0)
- # fid.close()
return fid, tree, directory
diff --git a/fiff/tag.py b/fiff/tag.py
index 9655cf3..d932272 100644
--- a/fiff/tag.py
+++ b/fiff/tag.py
@@ -129,7 +129,6 @@ def read_tag(fid, pos=None):
# All other data types
# Simple types
-
if tag.type == FIFF.FIFFT_BYTE:
tag.data = np.fromfile(fid, dtype=">B1", count=tag.size)
elif tag.type == FIFF.FIFFT_SHORT:
diff --git a/fiff/tests/test_evoked.py b/fiff/tests/test_evoked.py
new file mode 100644
index 0000000..52e7178
--- /dev/null
+++ b/fiff/tests/test_evoked.py
@@ -0,0 +1,32 @@
+import os
+import os.path as op
+
+import numpy as np
+from numpy.testing import assert_array_almost_equal, assert_equal
+
+import fiff
+
+MNE_SAMPLE_DATASET_PATH = os.getenv('MNE_SAMPLE_DATASET_PATH')
+fname = op.join(MNE_SAMPLE_DATASET_PATH, 'MEG', 'sample',
+ 'sample_audvis-ave.fif')
+
+def test_io_cov():
+ """Test IO for noise covariance matrices
+ """
+ data = fiff.read_evoked(fname)
+
+ fiff.write_evoked('evoked.fif', data)
+ data2 = fiff.read_evoked('evoked.fif')
+
+ print assert_array_almost_equal(data['evoked']['epochs'],
+ data2['evoked']['epochs'])
+ print assert_array_almost_equal(data['evoked']['times'],
+ data2['evoked']['times'])
+ print assert_equal(data['evoked']['nave'],
+ data2['evoked']['nave'])
+ print assert_equal(data['evoked']['aspect_kind'],
+ data2['evoked']['aspect_kind'])
+ print assert_equal(data['evoked']['last'],
+ data2['evoked']['last'])
+ print assert_equal(data['evoked']['first'],
+ data2['evoked']['first'])
diff --git a/fiff/tree.py b/fiff/tree.py
index 2ce7871..86b23c1 100644
--- a/fiff/tree.py
+++ b/fiff/tree.py
@@ -73,22 +73,22 @@ def make_dir_tree(fid, directory, start=0, indent=0, verbose=False):
# Add the id information if available
if block == 0:
- if directory[this].kind == FIFF_FILE_ID:
- tag = read_tag(fid, directory[this].pos)
- tree.id = tag.data
+ if directory[this].kind == FIFF_FILE_ID:
+ tag = read_tag(fid, directory[this].pos)
+ tree.id = tag.data
else:
- if directory[this].kind == FIFF_BLOCK_ID:
- tag = read_tag(fid, directory[this].pos)
- tree.id = tag.data
- elif directory[this].kind == FIFF_PARENT_BLOCK_ID:
- tag = read_tag(fid, directory[this].pos)
- tree.parent_id = tag.data
+ if directory[this].kind == FIFF_BLOCK_ID:
+ tag = read_tag(fid, directory[this].pos)
+ tree.id = tag.data
+ elif directory[this].kind == FIFF_PARENT_BLOCK_ID:
+ tag = read_tag(fid, directory[this].pos)
+ tree.parent_id = tag.data
this += 1
# Eliminate the empty directory
if tree.nent == 0:
- tree.directory = None
+ tree.directory = None
if verbose:
print '\t'*(indent+1) + 'block = %d nent = %d nchild = %d' % (
@@ -97,3 +97,58 @@ def make_dir_tree(fid, directory, start=0, indent=0, verbose=False):
last = this
return tree, last
+
+###############################################################################
+# Writing
+
+import numpy as np
+import struct
+from .constants import FIFF
+from .tag import Tag
+from .write import write_id, start_block, end_block, _write
+
+
+def copy_tree(fidin, in_id, nodes, fidout):
+ """
+ %
+ % fiff_copy_tree(fidin, in_id, nodes, fidout)
+ %
+ % Copies directory subtrees from fidin to fidout
+ %
+ """
+
+ if len(nodes) <= 0:
+ return
+
+ if not isinstance(nodes, list):
+ nodes = [nodes]
+
+ for node in nodes:
+ start_block(fidout, node.block)
+ if node['id'] is not None:
+ if in_id is not None:
+ write_id(fidout, FIFF.FIFF_PARENT_FILE_ID, in_id)
+
+ write_id(fidout, FIFF.FIFF_BLOCK_ID)
+ write_id(fidout, FIFF.FIFF_PARENT_BLOCK_ID, node['id'])
+
+ for d in node.directory:
+ # Do not copy these tags
+ if d.kind == FIFF.FIFF_BLOCK_ID or \
+ d.kind == FIFF.FIFF_PARENT_BLOCK_ID or \
+ d.kind == FIFF.FIFF_PARENT_FILE_ID:
+ continue
+
+ # Read and write tags, pass data through transparently
+ fidin.seek(d.pos, 0)
+
+ s = fidin.read(4*4)
+ tag = Tag(*struct.unpack(">iIii", s))
+ tag.data = np.fromfile(fidin, dtype='>B', count=tag.size)
+
+ _write(fidout, tag.data, tag.kind, 1, tag.type, '>B')
+
+ for child in node['children']:
+ copy_tree(fidin, in_id, child, fidout)
+
+ end_block(fidout, node.block)
diff --git a/fiff/write.py b/fiff/write.py
index e7e55c9..a0091b8 100644
--- a/fiff/write.py
+++ b/fiff/write.py
@@ -1,5 +1,7 @@
import time
+import array
import numpy as np
+from scipy import linalg
from .constants import FIFF
@@ -8,6 +10,8 @@ def _write(fid, data, kind, data_size, FIFFT_TYPE, dtype):
FIFFV_NEXT_SEQ = 0
if isinstance(data, np.ndarray):
data_size *= data.size
+ if isinstance(data, str):
+ data_size *= len(data)
fid.write(np.array(kind, dtype='>i4').tostring())
fid.write(np.array(FIFFT_TYPE, dtype='>i4').tostring())
fid.write(np.array(data_size, dtype='>i4').tostring())
@@ -79,7 +83,7 @@ def write_string(fid, kind, data):
%
"""
FIFFT_STRING = 10
- data_size = len(data)
+ data_size = 1
_write(fid, data, kind, data_size, FIFFT_STRING, '>c')
@@ -98,33 +102,39 @@ def write_name_list(fid, kind, data):
write_string(fid, kind, ':'.join(data))
-def write_float_matrix(fid, kind, data):
+def write_float_matrix(fid, kind, mat):
"""
%
% fiff_write_float_matrix(fid,kind,mat)
- %
+ %
% Writes a single-precision floating-point matrix tag
%
% fid An open fif file descriptor
% kind The tag kind
- % data The data matrix
+ % mat The data matrix
%
"""
-
FIFFT_FLOAT = 4
FIFFT_MATRIX = 1 << 30
FIFFT_MATRIX_FLOAT = FIFFT_FLOAT | FIFFT_MATRIX
- data_size = 4*data.size + 4*3
+ FIFFV_NEXT_SEQ = 0
- _write(fid, data, kind, data_size, FIFFT_MATRIX_FLOAT, '>f4')
+ data_size = 4 * mat.size + 4 * 3
- dims = np.empty(3, dtype=np.int)
- dims[0] = data.shape[1]
- dims[1] = data.shape[0]
- dims[3] = 2
+ fid.write(np.array(kind, dtype='>i4').tostring())
+ fid.write(np.array(FIFFT_MATRIX_FLOAT, dtype='>i4').tostring())
+ fid.write(np.array(data_size, dtype='>i4').tostring())
+ fid.write(np.array(FIFFV_NEXT_SEQ, dtype='>i4').tostring())
+ fid.write(np.array(mat, dtype='>f4').tostring())
+
+ dims = np.empty(3, dtype=np.int32)
+ dims[0] = mat.shape[1]
+ dims[1] = mat.shape[0]
+ dims[2] = 2
fid.write(np.array(dims, dtype='>i4').tostring())
+
def write_id(fid, kind, id_=None):
"""
%
@@ -142,7 +152,7 @@ def write_id(fid, kind, id_=None):
if id_ is None:
id_ = dict()
- id_['version'] = (1 << 16) | 2 # Version (1 << 16) | 2
+ id_['version'] = (1 << 16) | 2
id_['machid'] = 65536 * np.random.rand(2) # Machine id (andom for now)
id_['secs'] = time.time()
id_['usecs'] = 0 # Do not know how we could get this XXX
@@ -233,11 +243,195 @@ def end_file(fid):
% fid An open fif file descriptor
%
"""
-
data_size = 0
-
fid.write(np.array(FIFF.FIFF_NOP, dtype='>i4').tostring())
fid.write(np.array(FIFF.FIFFT_VOID, dtype='>i4').tostring())
fid.write(np.array(data_size, dtype='>i4').tostring())
fid.write(np.array(FIFF.FIFFV_NEXT_NONE, dtype='>i4').tostring())
fid.close()
+
+
+def write_coord_trans(fid, trans):
+ """
+ #
+ # fiff_write_coord_trans(fid,trans)
+ #
+ # Writes a coordinate transformation structure
+ #
+ # fid An open fif file descriptor
+ # trans The coordinate transfomation structure
+ #
+ """
+
+ FIFF_COORD_TRANS = 222
+ FIFFT_COORD_TRANS_STRUCT = 35
+ FIFFV_NEXT_SEQ = 0
+
+ #?typedef struct _fiffCoordTransRec {
+ # fiff_int_t from; /*!< Source coordinate system. */
+ # fiff_int_t to; /*!< Destination coordinate system. */
+ # fiff_float_t rot[3][3]; /*!< The forward transform (rotation part) */
+ # fiff_float_t move[3]; /*!< The forward transform (translation part) */
+ # fiff_float_t invrot[3][3]; /*!< The inverse transform (rotation part) */
+ # fiff_float_t invmove[3]; /*!< The inverse transform (translation part) */
+ #} *fiffCoordTrans, fiffCoordTransRec; /*!< Coordinate transformation descriptor */
+
+ data_size = 4*2*12 + 4*2
+ fid.write(np.array(FIFF_COORD_TRANS, dtype='>i4').tostring())
+ fid.write(np.array(FIFFT_COORD_TRANS_STRUCT, dtype='>i4').tostring())
+ fid.write(np.array(data_size, dtype='>i4').tostring())
+ fid.write(np.array(FIFFV_NEXT_SEQ, dtype='>i4').tostring())
+ fid.write(np.array(trans['from_'], dtype='>i4').tostring())
+ fid.write(np.array(trans['to'], dtype='>i4').tostring())
+
+ # The transform...
+ rot = trans['trans'][:3, :3]
+ move = trans['trans'][:3, 3]
+ fid.write(np.array(rot, dtype='>f4').tostring())
+ fid.write(np.array(move, dtype='>f4').tostring())
+
+ # ...and its inverse
+ trans_inv = linalg.inv(trans.trans)
+ rot = trans_inv[:3, :3]
+ move = trans_inv[:3, 3]
+ fid.write(np.array(rot, dtype='>f4').tostring())
+ fid.write(np.array(move, dtype='>f4').tostring())
+
+
+def write_ch_info(fid, ch):
+ """
+ %
+ % fiff_write_ch_info(fid,ch)
+ %
+ % Writes a channel information record to a fif file
+ %
+ % fid An open fif file descriptor
+ % ch The channel information structure to write
+ %
+ % The type, cal, unit, and pos members are explained in Table 9.5
+ % of the MNE manual
+ %
+ """
+
+ FIFF_CH_INFO = 203
+ FIFFT_CH_INFO_STRUCT = 30
+ FIFFV_NEXT_SEQ = 0
+
+ #typedef struct _fiffChPosRec {
+ # fiff_int_t coil_type; /*!< What kind of coil. */
+ # fiff_float_t r0[3]; /*!< Coil coordinate system origin */
+ # fiff_float_t ex[3]; /*!< Coil coordinate system x-axis unit vector */
+ # fiff_float_t ey[3]; /*!< Coil coordinate system y-axis unit vector */
+ # fiff_float_t ez[3]; /*!< Coil coordinate system z-axis unit vector */
+ #} fiffChPosRec,*fiffChPos; /*!< Measurement channel position and coil type */
+
+
+ #typedef struct _fiffChInfoRec {
+ # fiff_int_t scanNo; /*!< Scanning order # */
+ # fiff_int_t logNo; /*!< Logical channel # */
+ # fiff_int_t kind; /*!< Kind of channel */
+ # fiff_float_t range; /*!< Voltmeter range (only applies to raw data ) */
+ # fiff_float_t cal; /*!< Calibration from volts to... */
+ # fiff_ch_pos_t chpos; /*!< Channel location */
+ # fiff_int_t unit; /*!< Unit of measurement */
+ # fiff_int_t unit_mul; /*!< Unit multiplier exponent */
+ # fiff_char_t ch_name[16]; /*!< Descriptive name for the channel */
+ #} fiffChInfoRec,*fiffChInfo; /*!< Description of one channel */
+
+ data_size = 4*13 + 4*7 + 16;
+
+ fid.write(np.array(FIFF_CH_INFO, dtype='>i4').tostring())
+ fid.write(np.array(FIFFT_CH_INFO_STRUCT, dtype='>i4').tostring())
+ fid.write(np.array(data_size, dtype='>i4').tostring())
+ fid.write(np.array(FIFFV_NEXT_SEQ, dtype='>i4').tostring())
+
+ # Start writing fiffChInfoRec
+ fid.write(np.array(ch['scanno'], dtype='>i4').tostring())
+ fid.write(np.array(ch['logno'], dtype='>i4').tostring())
+ fid.write(np.array(ch['kind'], dtype='>i4').tostring())
+ fid.write(np.array(ch['range'], dtype='>f4').tostring())
+ fid.write(np.array(ch['cal'], dtype='>f4').tostring())
+ fid.write(np.array(ch['coil_type'], dtype='>i4').tostring())
+ fid.write(np.array(ch['loc'], dtype='>f4').tostring()) # writing 12 values
+
+ # unit and unit multiplier
+ fid.write(np.array(ch['unit'], dtype='>i4').tostring())
+ fid.write(np.array(ch['unit_mul'], dtype='>i4').tostring())
+
+ # Finally channel name
+ if len(ch['ch_name']):
+ ch_name = ch['ch_name'][:15]
+ else:
+ ch_name = ch['ch_name']
+
+ fid.write(np.array(ch_name, dtype='>c').tostring())
+ if len(ch_name) < 16:
+ dum = array.array('c', '\0' * (16 - len(ch_name)))
+ dum.tofile(fid)
+
+
+def write_dig_point(fid, dig):
+ """
+ %
+ % fiff_write_dig_point(fid,dig)
+ %
+ % Writes a digitizer data point into a fif file
+ %
+ % fid An open fif file descriptor
+ % dig The point to write
+ %
+ """
+
+ FIFF_DIG_POINT = 213
+ FIFFT_DIG_POINT_STRUCT = 33
+ FIFFV_NEXT_SEQ = 0
+
+ #?typedef struct _fiffDigPointRec {
+ # fiff_int_t kind; /*!< FIFF_POINT_CARDINAL,
+ # * FIFF_POINT_HPI, or
+ # * FIFF_POINT_EEG */
+ # fiff_int_t ident; /*!< Number identifying this point */
+ # fiff_float_t r[3]; /*!< Point location */
+ #} *fiffDigPoint,fiffDigPointRec; /*!< Digitization point description */
+
+ data_size = 5 * 4
+
+ fid.write(np.array(FIFF_DIG_POINT, dtype='>i4').tostring())
+ fid.write(np.array(FIFFT_DIG_POINT_STRUCT, dtype='>i4').tostring())
+ fid.write(np.array(data_size, dtype='>i4').tostring())
+ fid.write(np.array(FIFFV_NEXT_SEQ, dtype='>i4').tostring())
+
+ # Start writing fiffDigPointRec
+ fid.write(np.array(dig['kind'], dtype='>i4').tostring())
+ fid.write(np.array(dig['ident'], dtype='>i4').tostring())
+ fid.write(np.array(dig['r'][:3], dtype='>f4').tostring())
+
+
+def write_named_matrix(fid, kind, mat):
+ """
+ %
+ % fiff_write_named_matrix(fid,kind,mat)
+ %
+ % Writes a named single-precision floating-point matrix
+ %
+ % fid An open fif file descriptor
+ % kind The tag kind to use for the data
+ % mat The data matrix
+ %
+ """
+
+ start_block(fid, FIFF.FIFFB_MNE_NAMED_MATRIX)
+ write_int(fid, FIFF.FIFF_MNE_NROW, mat['nrow'])
+ write_int(fid, FIFF.FIFF_MNE_NCOL, mat['ncol'])
+
+ if len(mat['row_names']) > 0:
+ write_name_list(fid, FIFF.FIFF_MNE_ROW_NAMES, mat['row_names'])
+
+ if len(mat['col_names']) > 0:
+ write_name_list(fid, FIFF.FIFF_MNE_COL_NAMES, mat['col_names'])
+
+ write_float_matrix(fid,kind, mat.data)
+ end_block(fid, FIFF.FIFFB_MNE_NAMED_MATRIX)
+
+ return;
+
--
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