[med-svn] [python-mne] 168/353: maxfilter, EOG SSP computation
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:51 UTC 2015
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 <mluessi at nmr.mgh.harvard.edu>
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
More information about the debian-med-commit
mailing list