This is an automated email from the git hooks/post-receive script. yoh pushed a commit to tag 0.4 in repository python-mne.
commit 6c851019cc8b33146f596a985a4af61c4a62a9bf Author: Martin Luessi <[email protected]> Date: Mon Apr 30 14:52:11 2012 -0400 maxfilter, EOG SSP computation --- bin/mne_compute_proj_eog.py | 128 +++++++++++++++++++++++++++++++ mne/preprocessing/__init__.py | 3 +- mne/preprocessing/maxfilter.py | 170 +++++++++++------------------------------ mne/preprocessing/ssp.py | 82 +++++++++++++++++++- 4 files changed, 257 insertions(+), 126 deletions(-) diff --git a/bin/mne_compute_proj_eog.py b/bin/mne_compute_proj_eog.py new file mode 100755 index 0000000..cfaf345 --- /dev/null +++ b/bin/mne_compute_proj_eog.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python +"""Compute SSP/PCA projections for EOG artifacts + +You can do for example: + +$mne_compute_proj_eog.py -i sample_audvis_raw.fif --l-freq 1 --h-freq 100 --rej-grad 3000 --rej-mag 4000 --rej-eeg 100 +""" + +# Authors : Alexandre Gramfort, Ph.D. +# Martin Luessi, Ph.D. + +import sys +import os +import mne + + +if __name__ == '__main__': + + from optparse import OptionParser + + parser = OptionParser() + parser.add_option("-i", "--in", dest="raw_in", + help="Input raw FIF file", metavar="FILE") + parser.add_option("--tmin", dest="tmin", + help="Time before event in seconds", + default=-0.15) + parser.add_option("--tmax", dest="tmax", + help="Time after event in seconds", + default=0.15) + parser.add_option("-g", "--n-grad", dest="n_grad", + help="Number of SSP vectors for gradiometers", + default=2) + parser.add_option("-m", "--n-mag", dest="n_mag", + help="Number of SSP vectors for magnetometers", + default=2) + parser.add_option("-e", "--n-eeg", dest="n_eeg", + help="Number of SSP vectors for EEG", + default=2) + parser.add_option("--l-freq", dest="l_freq", + help="Filter low cut-off frequency in Hz", + default=5) + parser.add_option("--h-freq", dest="h_freq", + help="Filter high cut-off frequency in Hz", + default=35) + parser.add_option("-p", "--preload", dest="preload", + help="Temporary file used during computaion", + default='tmp.mmap') + parser.add_option("-a", "--average", dest="average", action="store_true", + help="Compute SSP after averaging", + default=False) + parser.add_option("--filtersize", dest="filter_length", + help="Number of taps to use for filtering", + default=2048) + parser.add_option("-j", "--n-jobs", dest="n_jobs", + help="Number of jobs to run in parallel", + default=1) + parser.add_option("--rej-grad", dest="rej_grad", + help="Gradiometers rejection parameter in fT/cm (peak to peak amplitude)", + default=2000) + parser.add_option("--rej-mag", dest="rej_mag", + help="Magnetometers rejection parameter in fT (peak to peak amplitude)", + default=3000) + parser.add_option("--rej-eeg", dest="rej_eeg", + help="EEG rejection parameter in uV (peak to peak amplitude)", + default=50) + parser.add_option("--rej-eog", dest="rej_eog", + help="EOG rejection parameter in uV (peak to peak amplitude)", + default=250) + parser.add_option("--avg-ref", dest="avg_ref", action="store_true", + help="Add EEG average reference proj", + default=False) + parser.add_option("--existing", dest="include_existing", action="store_true", + help="Inlucde the SSP projectors currently in the fiff file", + default=True) + parser.add_option("--bad", dest="bad_fname", + help="Text file containing bad channels list (one per line)", + default=None) + + options, args = parser.parse_args() + + raw_in = options.raw_in + + if raw_in is None: + parser.print_help() + sys.exit(-1) + + tmin = options.tmin + tmax = options.tmax + n_grad = options.n_grad + n_mag = options.n_mag + n_eeg = options.n_eeg + l_freq = options.l_freq + h_freq = options.h_freq + average = options.average + preload = options.preload + filter_length = options.filter_length + n_jobs = options.n_jobs + reject = dict(grad=1e-13 * float(options.rej_grad), + mag=1e-15 * float(options.rej_mag), + eeg=1e-6 * float(options.rej_eeg), + eog=1e-6 * float(options.rej_eog)) + avg_ref = options.avg_ref + include_existing = options.include_existing + bad_fname = options.bad_fname + + if bad_fname is not None: + bads = [w.rstrip().split()[0] for w in open(bad_fname).readlines()] + print 'Bad channels read : %s' % bads + else: + bads = [] + + if raw_in.endswith('_raw.fif') or raw_in.endswith('-raw.fif'): + prefix = raw_in[:-8] + else: + prefix = raw_in[:-4] + + eog_event_fname = prefix + '_eog-eve.fif' + + if average: + eog_proj_fname = prefix + '_eog_avg_proj.fif' + else: + eog_proj_fname = prefix + '_eog_proj.fif' + + mne.preprocessing.compute_proj_eog(raw_in, tmin, tmax, + n_grad, n_mag, n_eeg, l_freq, h_freq, average, preload, + filter_length, n_jobs, reject, bads, + avg_ref, include_existing, eog_proj_fname, eog_event_fname) + diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index 3dc95f2..41f62ea 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -6,4 +6,5 @@ # # License: BSD (3-clause) -from . ssp import compute_proj_ecg +from . maxfilter import apply_maxfilter +from . ssp import compute_proj_ecg, compute_proj_eog diff --git a/mne/preprocessing/maxfilter.py b/mne/preprocessing/maxfilter.py index 00a9d1c..35d86bd 100644 --- a/mne/preprocessing/maxfilter.py +++ b/mne/preprocessing/maxfilter.py @@ -4,7 +4,7 @@ # # License: BSD (3-clause) -import subprocess +import os from warnings import warn import numpy as np @@ -82,14 +82,11 @@ def fit_sphere_to_headshape(info): def apply_maxfilter(in_fname, out_fname, origin=None, frame='device', - in_order=None, out_order=None, bad=None, autobad=None, - badlimit=None, skip=None, data_format=None, force=False, - st=False, st_buflen=None, st_corr=None, mv_trans=None, + bad=None, autobad='off', skip=None, force=False, + st=False, st_buflen=16.0, st_corr=0.96, mv_trans=None, mv_comp=False, mv_headpos=False, mv_hp=None, - mv_hpistep=None, mv_hpisubt=None, mv_hpicons=False, - linefreq=None, lpfilt=None, site=None, cal=None, - ctc=None, magbad=False, regularize=None, iterate=None, - ds=None): + mv_hpistep=None, mv_hpisubt=None, mv_hpicons=True, + linefreq=None, mx_args='', overwrite=True): """ Apply NeuroMag MaxFilter to raw data. @@ -103,45 +100,33 @@ def apply_maxfilter(in_fname, out_fname, origin=None, frame='device', out_fname: string Output file name - origin: ndarray + origin: array-like or string Head origin in mm. If None it will be estimated from headshape points. frame: string ('device' or 'head') Coordinate frame for head center - in_order: int (or None) - Order of the inside expansion (None: use default) - - out_order: int (or None) - Order of the outside expansion (None: use default) - bad: string (or None) List of static bad channels (logical chnos, e.g.: 0323 1042 2631) - autobad: string ('on', 'off', 'n') (or None) - Sets automated bad channel detection on or off (None: use default) - - badlimit: float (or None) - Threshold for bad channel detection (>ave+x*SD) (None: use default) + autobad: string ('on', 'off', 'n') + Sets automated bad channel detection on or off skip: string (or None) Skips raw data sequences, time intervals pairs in sec, e.g.: 0 30 120 150 - data_format: string ('short', 'long', 'float') (or None) - Output data format (None: use default) - force: bool Ignore program warnings st: bool Apply the time-domain MaxST extension - st_buflen: float (or None) - MaxSt buffer length in sec (None: use default) + st_buflen: float + MaxSt buffer length in sec (disabled if st is False) - st_corr: float (or None) - MaxSt subspace correlation limit (None: use default) + st_corr: float + MaxSt subspace correlation limit (disabled if st is False) mv_trans: string (filename or 'default') (or None) Transforms the data into the coil definitions of in_fname, or into the @@ -152,51 +137,32 @@ def apply_maxfilter(in_fname, out_fname, origin=None, frame='device', mv_headpos: bool Estimates and stores head position parameters, but does not compensate - movements + movements (disabled if mv_comp is False) mv_hp: string (or None) Stores head position data in an ascii file + (disabled if mv_comp is False) mv_hpistep: float (or None) - Sets head position update interval in ms (None: use default) + Sets head position update interval in ms (disabled if mv_comp is False) mv_hpisubt: string ('amp', 'base', 'off') (or None) Subtracts hpi signals: sine amplitudes, amp + baseline, or switch off - (None: use default) + (disabled if mv_comp is False) mv_hpicons: bool Check initial consistency isotrak vs hpifit + (disabled if mv_comp is False) linefreq: int (50, 60) (or None) Sets the basic line interference frequency (50 or 60 Hz) (None: do not use line filter) - lpfilt: float (or None or True) - Corner frequency for IIR low-pass filtering - (None: don't use option: True: use default frequency) - - site: string (or None) - Loads sss_cal_<name>.dat and ct_sparse_<name>.fif - - cal: string (filename or 'off') (or None) - Uses the fine-calibration in <filename>, or switch off. + mx_args: string + Additional command line arguments to pass to MaxFilter - ctc: string (filename or 'off') (or None) - Uses the cross-talk matrix in <filename>, or switch off - - magbad: bool - Marks all magnetometers bad - - regularize: bool (or None) - Sets the component selection on or off (None: use default) - - iterate: int (or None) - Uses iterative pseudo-inverse, n iteration loops default n=10; - n=0 forces direct pseudo-inverse. (None: use default) - - ds: int (or None) - Applies downsampling with low-pass FIR filtering f is optional - downsampling factor (None: don't use option) + overwrite: bool + Overwrite output file if it already exists Returns ------- @@ -207,22 +173,19 @@ def apply_maxfilter(in_fname, out_fname, origin=None, frame='device', # check for possible maxfilter bugs def _mxwarn(msg): warn('Possible MaxFilter bug: %s, more info: ' - 'http://imaging.mrc-cbu.cam.ac.uk/meg/maxbugs') + 'http://imaging.mrc-cbu.cam.ac.uk/meg/maxbugs' % msg) - if mv_trans is not None and mv_comp is not None: + if mv_trans is not None and mv_comp: _mxwarn("Don't use '-trans' with head-movement compensation " "'-movecomp'") - if autobad is not None and (mv_headpos is not None or mv_comp is not None): + if autobad != 'off' and (mv_headpos or mv_comp): _mxwarn("Don't use '-autobad' with head-position estimation " "'-headpos' or movement compensation '-movecomp'") - if st and autobad is not None: + if st and autobad != 'off': _mxwarn("Don't use '-autobad' with '-st' option") - if lpfilt is not None: - _mxwarn("Don't use '-lpfilt' to low-pass filter data") - # determine the head origin if necessary if origin is None: print 'Estimating head origin from headshape points..' @@ -237,100 +200,59 @@ def apply_maxfilter(in_fname, out_fname, origin=None, frame='device', else: RuntimeError('invalid frame for origin') - # format command - cmd = ('maxfilter -f %s -o %s -frame %s -origin %d %d %d ' - % (in_fname, out_fname, frame, origin[0], origin[1], origin[2])) - - if in_order is not None: - cmd += '-in %d ' % in_order + if not isinstance(origin, str): + origin = '%0.1f %0.1f %0.1f' % (origin[0], origin[1], origin[2]) - if out_order is not None: - cmd += '-out %d ' % out_order + # format command + cmd = ('maxfilter -f %s -o %s -frame %s -origin %s ' + % (in_fname, out_fname, frame, origin)) if bad is not None: cmd += '-bad %s ' % bad - if autobad is not None: - cmd += '-autobad %s ' % autobad - - if badlimit is not None: - cmd += '-badlimit %0.4f ' % badlimit + cmd += '-autobad %s ' % autobad if skip is not None: cmd += '-skip %s ' % skip - if data_format is not None: - cmd += '-format %s ' % data_format - if force: cmd += '-force ' if st: cmd += '-st ' - - if st_buflen is not None: - if not st: - raise RuntimeError('st_buflen cannot be used if st == False') cmd += ' %d ' % st_buflen - - if st_corr is not None: - if not st: - raise RuntimeError('st_corr cannot be used if st == False') cmd += '-corr %0.4f ' % st_corr if mv_trans is not None: cmd += '-trans %s ' % mv_trans - if mv_comp is not None: + if mv_comp: cmd += '-movecomp ' if mv_comp == 'inter': cmd += ' inter ' - if mv_headpos: - cmd += '-headpos ' + if mv_headpos: + cmd += '-headpos ' - if mv_hp: - cmd += '-hp %s' % mv_hp + if mv_hp is not None: + cmd += '-hp %s' % mv_hp - if mv_hpisubt is not None: - cmd += 'hpisubt %s ' % mv_hpisubt + if mv_hpisubt is not None: + cmd += 'hpisubt %s ' % mv_hpisubt - if mv_hpicons: - cmd += '-hpicons ' + if mv_hpicons: + cmd += '-hpicons ' if linefreq is not None: cmd += '-linefreq %d ' % linefreq - if lpfilt is not None: - cmd += '-lpfilt ' - if not isinstance(lpfilt, bool): - cmd += '%0.1f ' % lpfilt - - if site is not None: - cmd += '-site %s ' % site - - if cal is not None: - cmd += '-cal %s ' % cal - - if ctc is not None: - cmd += '-ctc %s ' % ctc - - if magbad: - cmd += '-magbad ' - - if regularize is not None: - if regularize: - cmd += '-regularize on ' - else: - cmd += '-regularize off ' - - if iterate is not None: - cmd += '-iterate %d' % iterate + cmd += mx_args - if ds is not None: - cmd += '-ds %d ' % ds + if overwrite and os.path.exists(out_fname): + os.remove(out_fname) print 'Running MaxFilter: %s ' % cmd - out = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True) - print out + st = os.system(cmd) + if st != 0: + raise RuntimeError('MaxFilter returned non-zero exit status %d' % st) print '[done]' diff --git a/mne/preprocessing/ssp.py b/mne/preprocessing/ssp.py index 34adc41..01c4ba3 100644 --- a/mne/preprocessing/ssp.py +++ b/mne/preprocessing/ssp.py @@ -244,9 +244,89 @@ def compute_proj_ecg(in_fif_fname, tmin=-0.2, tmax=0.4, reject, bads, avg_ref, include_existing, ecg_proj_fname, ecg_event_fname) - return projs, ecg_events +def compute_proj_eog(in_fif_fname, tmin=-0.15, tmax=0.15, + n_grad=2, n_mag=2, n_eeg=2, l_freq=1.0, h_freq=35.0, + average=False, preload="tmp.mmap", + filter_length=2048, n_jobs=1, + reject=dict(grad=2000e-13, mag=3000e-15, eeg=50e-6, + eog=250e-6), bads=None, + avg_ref=False, include_existing=False, + ecg_proj_fname=None, ecg_event_fname=None): + """Compute SSP/PCA projections for EOG artifacts + + Parameters + ---------- + in_fif_fname: string + Input Raw FIF file + + tmin: float + Time before event in second + + tmax: float + Time after event in seconds + + n_grad: int + Number of SSP vectors for gradiometers + + n_mag: int + Number of SSP vectors for magnetometers + + n_eeg: int + Number of SSP vectors for EEG + + l_freq: float + Filter low cut-off frequency in Hz + + h_freq: float + Filter high cut-off frequency in Hz + + average: bool + Compute SSP after averaging + + preload: string (or True) + Temporary file used during computaion + + filter_length: int + Number of taps to use for filtering + + n_jobs: int + Number of jobs to run in parallel + + reject: dict + Epoch rejection configuration (see Epochs) + + bads: list + List with (additional) bad channels + + avg_ref: bool + Add EEG average reference proj + + include_existing: bool + Inlucde the SSP projectors currently in the fiff file + + eog_proj_fname: string (or None) + Filename to use for projectors (not saved if None) + + eog_event_fname: string (or None) + Filename to use for events (not saved if None) + + Returns + ------- + proj : list + Computed SSP projectors + + eog_events : ndarray + Detected ECG events + """ + + projs, eog_events = _compute_exg_proj('EOG', in_fif_fname, tmin, tmax, + n_grad, n_mag, n_eeg, l_freq, h_freq, + average, preload, filter_length, n_jobs, None, + reject, bads, avg_ref, include_existing, + ecg_proj_fname, ecg_event_fname) + return projs, eog_events -- 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
