Requires stormdb v.0.6.1 or greater (run on 0.7.dev0)
In [51]:
from os.path import join
proj_name = 'MEG_EEG-Training'
scratch_folder = join('/projects', proj_name, 'scratch')
mf_folder = join(scratch_folder, 'maxfilter', 'VS-meanpos') # for maxfilter output
scripts_folder = join('/projects', proj_name, 'scripts')
misc_folder = join('/projects', proj_name, 'misc')
trans_folder = join(scratch_folder, 'trans') # for transforms
In [2]:
from stormdb.access import Query
from stormdb.process import Maxfilter
from mne.io import Raw
from mne.bem import fit_sphere_to_headshape
from mne.transforms import rotation_angles, rotation3d, write_trans
import warnings
import os
import re
import numpy as np
In [3]:
# silence mne a bit
from mne.utils import set_log_level
set_log_level('ERROR')
In [44]:
tsss_buffer_len = 60.
tsss_corr_lim = 0.96 # don't expect this to make huge diff
# if you know that some channels are bad or flat, enter them here
# in the form ['2511', '2241']
mx_cmd = join(misc_folder, 'bin', 'maxfilter-2.2.15')
cal_db = join(misc_folder, 'databases', 'sss_cal.dat')
ctc_db = join(misc_folder, 'databases', 'ct_sparse.fif')
In [4]:
qr = Query(proj_name)
cur_sub = '0012'
Here we want to calculate the average initial head position and use movecomp to correct head motion to that origin.
In [15]:
description = 'run_*'
DATAblocks = qr.filter_series(description=description, subjects=cur_sub, modalities='MEG')
if len(DATAblocks) != 2:
raise RuntimeError('Not all 2 blocks found for {0}, please check!'.format(cur_sub))
for ib in range(len(DATAblocks)):
print('{:d}: {:s}'.format(ib + 1, DATAblocks[ib]['path']))
Get device to head-transformation from the initial HPI fit at the beginning of each acquisition. This consists of a translation and rotation, which are combined in to a single transformation matrix. Since both operations are affine transformations, we may simply average the initial matrices to obtain the mean head position and rotation.
In [39]:
init_xfm = []
init_rot = []
data_len = []
for bl in DATAblocks:
fname = join(bl['path'], bl['files'][0]) # first file is enough
with warnings.catch_warnings(): # suppress some annoying warnings for now
warnings.simplefilter("ignore")
raw = Raw(fname, preload=False, verbose=False)
data_len += [len(raw)/raw.info['sfreq']]
init_xfm += [raw.info['dev_head_t']['trans']]
# translations: info['dev_head_t']['trans'][:, 3][:-1]
init_rot += [raw.info['dev_head_t']['trans'][:3, :3]]
Find the average head position and calculate how far each block is from it (look for outliers).
In [24]:
mean_init_xfm = np.mean(np.stack(init_xfm), axis=0) # stack, then average over new dim
init_rot_angles = [rotation_angles(m) for m in init_rot]
mean_init_rot_xfm = rotation3d(*tuple(np.mean(np.stack(init_rot_angles),
axis=0))) # stack, then average, then make new xfm
assert(np.sum(mean_init_xfm[-1]) == 1.0) # sanity check result
mean_trans = info['dev_head_t'] # use the last info as a template
mean_trans['trans'] = mean_init_xfm # replace the transformation
mean_trans['trans'][:3, :3] = mean_init_rot_xfm # replace the rotation part
Plot the mean position in mm and each block's discrepancy from the mean.
In [25]:
mean_init_headpos = mean_trans['trans'][:-1, -1] # meters
print('Mean head position (device coords): ({:.1f}, {:.1f}, {:.1f}) mm'.\
format(*tuple(mean_init_headpos*1e3)))
print('Block discrepancies from mean:')
for ib, xfm in enumerate(init_xfm):
diff = 1e3 * (xfm[:-1, -1] - mean_init_headpos)
rmsdiff = np.linalg.norm(diff)
print('\tblock {:d}: norm {:.1f} mm ({:.1f}, {:.1f}, {:.1f}) mm '.\
format(ib + 1, rmsdiff, *tuple(diff)))
Do the same for rotations
In [26]:
mean_rots = rotation_angles(mean_trans['trans'][:3, :3]) # these are in radians
mean_rots_deg = tuple([180. * rot / np.pi for rot in mean_rots]) # convert to deg
print('Mean head rotations (around x, y & z axes): ({:.1f}, {:.1f}, {:.1f}) deg'.\
format(*mean_rots_deg))
print('Block discrepancies from mean:')
for ib, rot in enumerate(init_rot):
cur_rots = rotation_angles(rot)
diff = tuple([180. * cr / np.pi - mr for cr, mr in zip(cur_rots, mean_rots_deg)])
print('\tblock {:d}: ({:.1f}, {:.1f}, {:.1f}) deg '.\
format(ib + 1, *tuple(diff)))
Save the new mean transformation to a file to be used later in the maxfilter-option trans.
In [27]:
mean_trans_folder = join(trans_folder, cur_sub)
if not os.path.exists(mean_trans_folder):
os.makedirs(mean_trans_folder)
mean_trans_file = join(mean_trans_folder, 'mean-trans.fif')
write_trans(mean_trans_file, mean_trans)
In [28]:
set_log_level('INFO')
rad, origin_head, ori_dev = fit_sphere_to_headshape(info,
dig_kinds='extra',
units='mm')
set_log_level('ERROR')
In [68]:
mf = Maxfilter(proj_name, bad=sub_specific_bad_chans)
In [69]:
mfopts = dict(
origin = '{:.1f} {:.1f} {:.1f}'.format(*tuple(origin_head)), # mm
frame = 'head',
force = True, # overwrite if needed
autobad = 'on', # or use xscan first
st = True, # use tSSS
st_buflen = tsss_buffer_len, # parameter set in beg. of notebook
st_corr = tsss_corr_lim, # parameter set in beg. of notebook
movecomp = True,
trans = mean_trans_file, # compensate to mean initial head position (saved to file),
# use None for initial head position
logfile = None, # we replace this in each loop
hp = None, # head positions, replace in each loop
n_threads = 4 # antal kerner på isis, max 12, være solidarisk!
)
mne-python likes raw and raw-like (tsss) files that are part of a long (>2GB) continuous acquisition to follow the convention:
filename_raw_tsss.fif (first file)filename_raw_tsss-1.fif (second file, etc.)
In [70]:
out_folder = join(mf_folder,
'tsss_st{:.0f}_corr{:.0f}'.format(mfopts['st_buflen'],
np.round(100 * mfopts['st_corr'])),
cur_sub)
# Check that output path exists
if not os.path.exists(out_folder):
os.makedirs(out_folder)
In [71]:
for blockno, bl in enumerate(DATAblocks):
for fileno, fil in enumerate(bl['files']):
in_fname = join(bl['path'], bl['files'][fileno])
series_name = re.search('(run_1|run_2|run_3)',
bl['seriename']).group(1)
out_fname = join(out_folder, '{0}_raw_tsss.fif'.format(series_name))
if fileno > 0:
out_fname = out_fname[:-4] + '-{:d}.fif'.format(fileno)
mfopts['logfile'] = out_fname[:-3] + 'log'
mfopts['hp'] = out_fname[:-3] + 'pos'
mf.build_cmd(in_fname, out_fname, **mfopts)
In [72]:
mf.print_input_output_mapping()
If in doubt, uncomment this line to see the actual commands that will execute:
In [73]:
mf.commands
Out[73]:
In [61]:
# mf.cluster.get_load()
In [62]:
mf.submit()
In [66]:
mf.status
To kill a submitted (or even running) job
In [63]:
# mf.kill(3511711)
To kill all submitted jobs:
In [65]:
# mf.kill()
To kill a job in the shell:
qdel JOB_NUMBERor for all jobs (in your name):
qdel *
In [59]:
try:
mf_list
except NameError:
mf_list = []
mf_list += [mf]
In [60]:
for cur_mf in mf_list:
cur_mf.status
In [ ]: