#!/usr/bin/env python ID = "$Id: preprocess.py 198 2009-09-17 23:37:21Z jmo $"[1:-1] # Written by John Ollinger # # University of Wisconsin, 8/16/09 #Copyright (c) 2006-2007, John Ollinger, University of Wisconsin #All rights reserved. # #Redistribution and use in source and binary forms, with or without #modification, are permitted provided that the following conditions are #met: # # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # # * Redistributions in binary form must reproduce the above # copyright notice, this list of conditions and the following # disclaimer in the documentation and/or other materials provided # with the distribution. # #THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS #"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT #LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR #A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT #OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, #SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT #LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, #DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY #THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT #(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE #OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # # ** This software was designed to be used only for research purposes. ** # ** Clinical uses are not recommended, and have never been evaluated. ** # ** This software comes with no warranties of any kind whatsoever, ** # ** and may not be useful for anything. Use it at your own risk! ** # ** If these terms are not acceptable, you aren't allowed to use the code.** import sys import os from os import F_OK, R_OK, W_OK from optparse import OptionParser from file_io import Wimage, Header, exec_cmd, ispfile, isIfile, \ abspath_to_relpath, append_history_note from dicom import isdicom from numpy import zeros from datetime import datetime from wimage_lib import execCmd, ExecError, except_msg, \ send_email, get_tmp_space if 'linux' in sys.platform: import dl import pickle import yaml from yaml.scanner import ScannerError import smtplib from email.mime.text import MIMEText import socket from socket import gethostname import time import constants as c import traceback from stat import S_IRWXU, S_IRWXG from bz2 import BZ2File ID = "$Id: preprocess.py 198 2009-09-17 23:37:21Z jmo $"[1:-1] def echo_ID(): return ID #EXAMS_FILE='/study/scanner/preprocessed_exams.txt' DEFAULT_SKIP=5 ERROR=1 IGNORE=1 OK=0 UMASK_FILE=0113 UMASK_DIR=0002 ZIP_EXCLUDE = \ '-not \( -name "*.txt" \) ' + \ '-not \( -name "*.yaml" \) ' + \ '-not \( -name "*.doc" \) ' + \ '-not \( -name "*.xls" \) ' + \ '-not \( -name "*.ppt" \) ' + \ '-not \( -name "*.bsh" \) ' + \ '-not \( -name "*.jpg" \) ' + \ '-not \( -name "*.png" \) ' + \ '-not \( -name "*.pro" \) ' + \ '-not \( -name "*.es" \) ' + \ '-not \( -name "*.ebs" \) ' + \ '-not \( -name "*.pdf" \) ' + \ '-not \( -name "*.py" \) ' + \ '-not \( -name "*.sta" \) ' + \ '-not \( -name "*.m" \) ' + \ '-not \( -name "*.exe" \) ' + \ '-not \( -name "*.spo" \) ' + \ '-not \( -name "*.sav" \) ' + \ '-not \( -name "*.1D" \) ' + \ '-not \( -name "*.xcl" \) ' + \ '-not \( -name "*.pl" \) ' decompress_cmds = {'.gz':'gunzip --to-stdout', '.bz2':'bunzip2 --stdout'} #********************************************************************** def EmailResults(recipient, error_mesg, topdir, dumpfile, logfile): #********************************************************************** """ Email summary of results to user. """ if recipient is None: return elif 'noname' in recipient: return sender = 'preprocess' if 'Abnormal' in error_mesg > 0: subject = 'Problem while preprocessing %s' % topdir else: subject = 'Preprocessing complete for %s' % topdir mssg = error_mesg if logfile is not None and isinstance(logfile, str): f = open(logfile, 'r') lines = f.readlines() f.close() logged_errors = '' for i in xrange(len(lines)): if 'rror' in lines[i]: mssg += ''.join(lines[i-1:]) break if dumpfile is not None: f = open(dumpfile,'r') mssg += '\nSummary of processing:\n' mssg += f.read() f.close() send_email(recipient, subject, mssg, sender) #*********************** def hms_to_secs(hms_in): #*********************** # Converts time in the string format of hh:mm:ss to seconds since midnight. hms = hms_in.split(':') return float(hms[2]) + 60.*(float(hms[1]) + 60.*float(hms[0])) class FakeOptions(): def __init__(self): self.debug_tmp = False self.outdir = None self.align_epis = False self.skull_strip = False self.keep_epi_raw = False self.keep_epi_mot = False self.fake_opts = True self.epi_keys = {} self.exclude_paths = [] class DataInfo(object): def __init__(self,topdir, redo=False, verbose=False, dry_run=False, \ skip=None, scratch='/scratch', no_email=False): # Get the version of the software. try: f = open(c.FMRI_TOOLS_VERSION_FILE,'r') self.version = f.read().strip() f.close() except IOError: self.version = '0.0.0' self.scratch = scratch self.no_email = no_email self.prog = os.path.basename(sys.argv[0])[:-3].strip() self.verstring = "\n%s: version %s\n\n" % (self.prog, self.version) if self.prog == 'preprocess': print self.verstring self.anatref_entry = None self.ntype = {} self.tmplt = None self.info = {} self.fmaps = {} self.pfiles = [] self.epirt_paths = [] self.refdats = {} self.epi_series = [] self.pfiles_recon = [] self.topdir = os.path.abspath(topdir) self.logdir = None self.verbose = verbose self.anatomical = None self.anat_ref = None self.n_epi = 1 self.epi_info = {} self.found_data = False self.visited_dirs = set() self.errors = False # Get options. if not hasattr(self, 'opts'): self.opts = FakeOptions() else: self.opts.fake_opts = False # Create temporary directory. if self.opts.debug_tmp: self.tmpdir = '/tmp/debug_tmp' self.MakeDir(self.tmpdir) else: self.tmpdir = get_tmp_space(1000, c.backup_tmp_dir) # Check paths to be excluded. self.exclude_paths = [] for path in self.opts.exclude_paths: apath = os.path.abspath(path) if os.path.exists(apath): self.exclude_paths.append(os.path.abspath(path)) else: sys.stderr.write('\n*** Nonexistent path specified in --exclude argument.\n***Invalid path: %s\n***Aborting ...\n' % apath) self.term_mesg = '\nAbnormal termination.\n' self.CleanUp() self._ProcessTemplate(self.topdir) # Check for keywords defining which epis to reconstruct. self.epi_keys = {} for keyword in self.opts.epi_keys: wds = keyword.split(':') self.epi_keys[wds[0]] = wds[1] # Table mapping image types to entries in the master info table. self.entry_map = {'epi':[], 'anat':[], 'fmap':[], 'dti':[], 'first_epi':[]} self.type = {'efgre3d':'T1High', \ '3dir':'T1High', \ 'bravo':'T1High', \ 'fse-xl':'T2', \ 'frfseopt':'T2', \ 'cubet2':'T2', \ '2dfast':'fmap', \ 'epibold':'epi', \ 'fcmemp':'T1se', \ 'fse':'T2', \ '3-plane':'localizer', \ 'epi':'epi', \ 'epirt':'epi', \ 'epi2':'dti', \ 'dti':'dti'} self.GetInfoMethods = {\ 'T1High': self._AnatInfo, \ 'T2': self._AnatInfo, \ 'fse': self._AnatInfo, \ 'T1se': self._AnatInfo, \ 'fmap': self._FmapInfo, \ 'dti': self._DtiInfo, \ 'epi': self._EpiInfo} self.anat_types = ['T1High', 'T2', 'T1se', 'fse'] self.suffix = {\ 'brik':'+orig.BRIK', \ 'nii':'.nii', \ 'n+1':'.nii', \ 'ni1':'.hdr', \ 'ana':'.hdr'} self.infile_suffix = {\ 'brik':'+orig', \ 'nii':'.nii', \ 'n+1':'nii'} while self.topdir.endswith('/'): self.topdir = self.topdir[:-1] if hasattr(self, 'LogProcess'): self.LogProcess() # Look for data to process. self.WalkPath(self.topdir) if os.path.islink('%s/anatomicals' % self.topdir): # os.walk won't follow links, so do this one manually. pathname = os.readlink('%s/anatomicals' % self.topdir) self.WalkPath(pathname) # Pair-up fieldmaps with EPI's self._GetFmapInfo() # Assocate a ref.dat file with each EPI. self._GetRefdat() # Order the EPIs so the names are correct. self._GetEpiOrder() # Associate each EPI with an anatomical, determine if it was # acquired before or after the epi self._FindBaseEpi() self.DetermineBaseFrame() def Gunzip(self, pathname): cmd = "gunzip -rfq %s" % pathname try: execCmd(cmd) except ExecError: self.errors = True errstr = "%s\nCould not unzip data files in %s" % (cmd, pathname) +\ "\nChange permissions and rerun the script." raise ExecError(errstr) def Gzip(self, pathname): cmd = 'find %s -type f %s' % (pathname, ZIP_EXCLUDE) try: print 'Gzip starting time: %s' % time.strftime('%H:%M:%S') execCmd(cmd) print 'Gzip ending time: %s' % time.strftime('%H:%M:%S') except ExecError: self.errors = True print 'Gzip ending time: %s' % time.strftime('%H:%M:%S') errstr = "Could not gzip data files in %s" % pathname + \ "Change permissions and compress by hand." # Don't call this a fatal error. sys.stderr.write('%s\n' % errstr) raise ExecError(errstr) def CheckCompression(self, filename): if filename.endswith('.gz'): fname = filename[:-3] compression = '.gz' elif filename.endswith('.bz2'): fname = filename[:-4] compression = '.bz2' else: fname = filename compression = None return (fname, compression) def WalkPath(self, pathname): "Check subdirectories of topdir for data we can process.""" pfile_basenames = [] self.n_t1low = 1 self.n_t1high = 1 followed_links = {} for root, dnames, fnames in os.walk(pathname, topdown=True): # Ensure we never parse the same directory twice idx = 0 for dname in dnames: full_dname = '%s/%s' % (root, dname) if os.path.islink(full_dname): # Edit directories so links will be followed (once). linked_path = os.readlink(full_dname) if followed_links.has_key(linked_path): sys.stderr.write(\ 'Ignoring link to avoid infinite loop: %s\n' % dname) else: followed_links[linked_path] = True dnames[idx] = linked_path idx += 1 real_root = os.path.realpath(root) if real_root in self.visited_dirs: continue self.visited_dirs.add(real_root) self.current_entry = root if self.current_entry.endswith('.yaml'): self.current_entry = os.path.dirname(self.current_entry) if root.startswith(pathname + '/orig'): continue if self.verbose: print 'Checking %s ...' % root if len(fnames) > 1: dirname, filename = self._yaml_filename(root) # os.path.dirname(fnames[0])) fullpath = '%s/%s' % (dirname, filename) if os.path.exists(fullpath): fnames = [filename] + fnames for fname in fnames: if fname.endswith('.yaml'): # Move yaml file to the beginning of the list. fnames.remove(fname) fnames = [fname] + fnames for name in fnames: fname, compression = self.CheckCompression(name) fullpath = '%s/%s' % (root,name) # Check if file is in an excluded directory. exclude = False for expath in self.exclude_paths: if fullpath.startswith(expath): exclude = True sys.stdout.write('*** Excluding data: %s\n' % \ fullpath) break if exclude: break if not os.path.isdir(fullpath) and os.path.exists(fullpath) and\ not os.path.basename(fullpath) in pfile_basenames: # Interrogate header of file specified by fullpath. info = self._GetImageInfo(fullpath) if info != None: if info['type'] == 'refdat': continue elif info['type'] == 'break': break info['compression'] = compression info['pfile_decomp'] = '%s/%s' % (root, fname) self.found_data = True if info['type'] == 'localizer': break elif info['data_filetype'] == 'ge_data': self.info[fullpath] = info pfile_basenames.append(os.path.basename(fullpath)) elif os.path.isdir(root): self.info[root] = info else: self.info[os.path.dirname(root)] = info if info['data_filetype'] == 'dicom' or \ info['data_filetype'] == 'ge_ifile': # Only look at one header per series or run. break # Prune duplicate p-files. pfiles = {} for fname in self.pfiles: # Non-unique entries will be overwritten. pfiles[os.path.basename(fname)] = fname self.pfiles = pfiles.values() def _FindTemplateFile(self, topdir): fnames = os.listdir(topdir) for fname in fnames: filename = '%s/%s' % (topdir, fname) if filename.endswith('.yaml') and not os.path.isdir(filename) and \ os.path.exists(filename): f = open(filename, 'r') magic_code = f.read(22) f.close() if '#!fmri_file_template' in magic_code: return filename return None def _LoadTemplate(self,fname): f = open(fname, 'r') data = f.read() data = data.replace('\t',' ') if '\t' in data: errstr = \ 'Illegal tabs encountered in template file. Use spaces instead.' raise ScannerError(errstr) proc.LogErrors(errstr) tmplt = yaml.load(data) f.close() return tmplt def _GetTemplate(self): # First read default template. tmplt = self._LoadTemplate(c.preproc_template_default) tmplt['proc'] = self.topdir self.template_type = 'default' if self.template_file is not None: tmplt.update(self._LoadTemplate(self.template_file)) self.template_type = 'command-line' found_template = True else: # Find a study specific template file. study_template_file = self._FindTemplateFile('%s/..' % self.topdir) if study_template_file is not None: # Merge study template into default, study template has precedence. if self.verbose: print "Using study template at " + study_template_file tmplt.update(self._LoadTemplate(study_template_file)) self.template_type = 'study-specific' found_template = True else: found_template = False # Now look for a subject-specific template file. subject_template_file = self._FindTemplateFile('%s' % self.topdir) if subject_template_file is not None: # Merge subject template, subject template has precedence. if self.verbose: print "Using subject-specific template at %s" % \ subject_template_file tmplt.update(self._LoadTemplate(subject_template_file)) self.template_type = 'study-specific' found_template = True if not found_template: raise RuntimeError('Could not find template file.') if tmplt.get('subject','same') == 'same': # Default subdirectory is same as data directory. tmplt['subject'] = self.topdir.split('/')[-1] else: if not isinstance(tmplt['subject'],str): errstr = 'preprocess: Invalid subject number. Be sure to ' + \ 'enclose the subject number item with double quotes.' raise RuntimeError(errstr) # Keys that apply to all EPIs. self.fsl_flip = tmplt.get('fsl_flip', False) if self.fsl_flip: self.flip_opts = '-LT' else: self.flip_opts = '' return tmplt def _ProcessTemplate(self,topdir): """ Set directory names. """ self.dicomdir = "%s/anatomicals" % self.topdir self.rawdir = "%s/raw" % topdir self.rawdirs = {} tmplt = self._GetTemplate() if self.opts.outdir is not None: # Override template output directory. tmplt['top_outdir'] = self.opts.outdir self.tmplt = tmplt if len(tmplt['top_outdir']) == 0: tmplt['top_outdir'] = self.topdir raise RuntimeError('Template file must specify an output directory.') # tmplt['subject'] = 'orig' self.procdir = os.path.abspath("%s/%s" % \ (tmplt['top_outdir'],tmplt['subject'])) self.anatdir = "%s/anat" % self.procdir self.anat_ref = tmplt['anat'].get('anat_ref', None) self.fmapdir = "%s/%s" % (self.procdir,tmplt['fmap']['outdir']) self.dtidir = "%s/%s" % (self.procdir,tmplt['dti']['outdir']) self.logdir = "%s/%s" % (self.procdir,tmplt['logdir']) self.skip = tmplt.get('skip', DEFAULT_SKIP) self.episetup_dir = "%s/%s" % (self.procdir,tmplt['first_epi']) self.fsl_cmpblty = tmplt.get('fsl_compatibility',False) self.epi_file_format = self.tmplt['epi_file_format'] self.censor_thresh = tmplt.get('censor_threshold', 2.) self.censor_interleave = tmplt.get('censor_interleave', True) # self.server_userid = self.tmplt.get('server_userid','default') # Overide flags for aligning EPIs and skull-stripping with command- # line options. if self.opts.align_epis: self.epi_align = True else: self.epi_align = self.tmplt.get('epi_align', False) if self.opts.skull_strip: self.skull_strip = True else: self.skull_strip = self.tmplt.get('skull_strip', False) # Create log file now so it can be used immediately. if not os.path.exists(self.logdir): if self.verbose: print 'mkdir %s' % self.logdir if not self.opts.fake_opts: self.MakeDir(self.logdir) self._ProcessTemplateEpiInfo() def _CheckForEmbeddedBlanks(string): if ' ' in string[1:-1]: # Embedded blank in filename. Missing comma? errstr = 'preprocess: Embedded blank in file name. Be sure '+\ 'that each file in the list is delimited by a comma.\n' raise RuntimeError(errstr) # def _ProcessTemplateEpiInfo(self): # Extract EPI groups and setup data structures used for sorting. self.epidirs = [] epi_info = {} keys = self.tmplt.keys() # Move epidir_dflt key to the end. keys.remove('epidir_dflt') keys = keys + ['epidir_dflt'] nepi_keys = 0 acq_orders = [] for key in self.tmplt.keys(): item = self.tmplt[key] if key == 'epidir_dflt': # Make sure this is always ordered last. self.tmplt[key]['acq_order'] = 1000 if nepi_keys > 0: # Don't process if there are explicitly definied epis. continue # For each group of epis, create a list containing an index to the # next name to be used and the orientation. if isinstance(item, dict): if not item.get('type','') == 'epi': continue nepi_keys += 1 epidir = '%s/%s' % (self.procdir, item['outdir']) if epidir not in self.epidirs: self.epidirs.append(epidir) plane = item.get('plane','any') acq_order = item.get('acq_order',0) # if epi_info.has_key(acq_order) and key != 'epidir_dflt': if acq_order in acq_orders: if key != 'epidir_dflt': raise RuntimeError(self.template_type + \ ' template has non-unique values of acq_order') elif epi_info.has_key(acq_order): # Overwrite lower-priority template file. del epi_info[acq_order] acq_orders.append(acq_order) epi_info[acq_order] = {'tmplt_key':key, \ 'plane':plane, \ 'anat_ref':[], \ 'subdir':[], \ 'names':[]} for name in item['names']: epi_info[acq_order]['names'].append(\ '%s/%s' % (epidir, name)) epi_info[acq_order]['anat_ref']. \ append(item.get('anat_ref','anat')) epi_info[acq_order]['subdir'].append(item['outdir']) epi_tmpdir = '%s/%s' % (self.tmpdir, item['outdir']) if epi_tmpdir not in self.epidirs: self.epidirs.append(epi_tmpdir) # Delete duplicate names, sort, and create final list. epi_acqs = epi_info.keys() epi_acqs.sort() self.epinames = {} for acq in epi_acqs: if len(epi_acqs) > 1 and \ epi_info[acq]['tmplt_key'] == 'epidir_dflt': # Delete epi group in default template if others were supplied. continue plane = epi_info[acq]['plane'] if not self.epinames.has_key(plane): self.epinames[plane] = \ {'n_epi':0, 'names':[], 'anat_ref':[],'subdir':[]} self.epinames[plane]['names'] += epi_info[acq]['names'] self.epinames[plane]['anat_ref'] += epi_info[acq]['anat_ref'] self.epinames[plane]['subdir'] += epi_info[acq]['subdir'] def _yaml_filename(self, path): # Synthesize yaml header filename from directory name. fullpath = os.path.abspath(path) if not os.path.isdir(fullpath): dirname = os.path.dirname(fullpath) else: dirname = path if dirname.endswith('/'): dirname = dirname[:-1] fname = dirname.split('/')[-1] + '.yaml' return dirname, fname def _FmapInfo(self, info, path): if not os.path.isdir(path): pathdir = os.path.dirname(path) else: pathdir = path fnames = os.listdir(pathdir) if len(fnames) < 8*self.hdr['zdim']: errstr = "\n*** Fieldmap data are incomplete: %s ***\n\n" % path self.LogErrors(errstr) return ERROR info['outdir'] = '%s/%s' % (self.procdir, self.tmplt['fmap']['outdir']) info['imgfile'] = '%s/fmap_%s' % (info['outdir'],info['plane'].strip()) info['magfile'] = '%s_mag' % info['imgfile'] info['filetype'] = 'nii' info['echo_spacing'] = self.tmplt['fmap']['echo_spacing'] if self.tmplt['fmap'].get('anat_ref','fieldmap') == 'fieldmap' or \ self.tmplt['fmap'].get('anat_ref','fieldmap') == 'default': # Default action is to align the fieldmap with the anatomical. This # attribute might not exist yet, so flag _GetFmapInfo to do after # all raw data has been surveyed. info['anat_ref'] = 'self.anatomical' else: info['anat_ref'] = self.tmplt['fmap']['anat_ref'] info['matfile'] = '%s_matfile.aff12.1D' % info['imgfile'] info['imgfile_m'] = info['imgfile'] self.fmaps[pathdir] = info['imgfile'] self.entry_map['fmap'].append(self.current_entry) return OK def _AnatInfo(self, info, path): """ Get T1 and T2 weighted structural image info.""" if info['data_filetype'] == 'ge_data': return ERROR info['outdir'] = '%s/%s' % (self.procdir, self.tmplt['anat']['outdir']) if info['psdname'] == 'efgre3d': if self.hdr['zsize'] > 1.25: # Only one slab acquired. Assume thick slices. if self.n_t1low < len(self.tmplt['anat'].get('t1_low',[])): name = self.tmplt['anat']['t1_low'][self.n_t1low-1] else: name = 'T1Low_%d' % self.n_t1low info['imgfile_m'] = '%s/%s_m' % (info['outdir'], name) self.n_t1low += 1 else: if not self.tmplt['anat'].has_key('t1_high'): self.tmplt['anat']['t1_high'] = ['T1High'] if self.n_t1high <= len(self.tmplt['anat']['t1_high']): name = self.tmplt['anat']['t1_high'][self.n_t1high-1] else: name = 'T1High_%d' % self.n_t1high info['imgfile_m'] = '%s/%s' % (info['outdir'], name) self.n_t1high += 1 if self.anat_ref == 'T1High': self.anat_ref = path else: psdname = info['psdname'] name = self.type.get(psdname, info['psdname']) if self.ntype.has_key(psdname): self.ntype[psdname] += 1 name = '%s_%0d' % (name, self.ntype[psdname]) else: self.ntype[psdname] = 1 if name == self.tmplt['anat'].get('anat_norm','T1High'): self.anatref_entry = self.current_entry info['filetype'] = self.tmplt['anat']['format'] info['imgfile'] = '%s/%s' % (info['outdir'], name) info['imgfile_skstrip'] = '%s_skstrip' % info['imgfile'] info['matfile'] = '%s_matfile.aff12.1D' % info['imgfile'] # name used by epi to refer to this anatomical. info['ref_name'] = name if self.tmplt.has_key('anat_norm'): # Use this image to normalize to atlas coordinates. info['anat_norm'] = '%s/%s' % \ (info['outdir'],self.tmplt['anat_norm']) else: info['anat_norm'] = '%s/T1High' % info['outdir'] self.anatomical = '%s%s' % \ (info['anat_norm'],self.suffix[info['filetype']]) # Save this entry in a table for later use. self.entry_map['anat'].append(self.current_entry) return OK def _DtiInfo(self, info, path): info['outdir'] = '%s/%s' %(self.procdir, self.tmplt['dti']['outdir']) info['filetype'] = self.tmplt['dti']['format'] info['pepolar'] = self.tmplt['dti']['pepolar'] info['imgfile'] = '%s/s%s_dti_%0ddir' % \ (info['outdir'],info['series'], info['tdim']-1) if not self.tmplt['dti'].has_key('anat_ref'): self.tmplt['dti']['anat_ref'] = 'fieldmap' if self.tmplt['dti']['anat_ref'] == 'fieldmap': # Default action is to align the fieldmap with the anatomical. info['anat_ref_entry'] = 'fmap' else: info['anat_ref'] = self.tmplt['dti']['anat_ref'] # Save this entry in a table for later use. self.entry_map['dti'].append(self.current_entry) return OK def _NoInfo(self, info, path): if info['type'] is None: return else: errstr = 'Invalid "type" specified.\ninfo: %s\n\n' % info['type'] self.LogErrors(errstr) raise RuntimeError(errstr) return OK def _EpiInfo(self, info, path): """ Create list of epis in pfile format (epi_series) and of epis in dicom format (epirt_paths) """ epi_vals = {'tdim':self.hdr['tdim'], 'plane':self.hdr['plane'], \ 'SeriesNumber':self.hdr['subhdr']['SeriesNumber']} for key in self.epi_keys.keys(): if self.epi_keys[key] != str(epi_vals[key]): # Return None, which will cause these data to be ignored. return None # Early versions of the EPIC software saved p-files for the setup # epis. Don't process these (or any epi with fewer than five frames). if self.hdr['tdim'] < 5 and info['data_filetype'] == 'ge_data': return None info['echo_spacing'] = self.shdr['EffEchoSpacing']/1000. if info['data_filetype'] == 'dicom': # Entry is name of dirctory for dicom images. if not os.path.isdir(path): entry = os.path.dirname(path) else: entry = path else: # Otherwise it is the name of a directory containing p-files. entry = path # if entry.endswith('.yaml'): or info['psdname'] == 'epirt': # entry = os.path.dirname(entry) if info['data_filetype'] == 'ge_data' and info['type'] is not None: # Found a pfile. Add it to the list. if entry not in self.pfiles and info['tdim'] > 2: self.pfiles.append(entry) self.entry_map['epi'].append(entry) if info['series'] not in self.epi_series: self.epi_series.append(info['series']) elif info['data_filetype'] == 'dicom' and \ info['psdname'] == 'epibold': # This is the initial EPI done during setup. info['outdir'] = self.episetup_dir info['type'] = 'first_epi' self.entry_map['first_epi'].append(entry) info['imgfile'] = '%s/first_epi_%d' % \ (self.episetup_dir, len(self.entry_map['first_epi'])) elif (info['psdname'] == 'epirt' or info['psdname'] == 'epi') \ and info['tdim'] > 2: # This is an epi reconstructed on the scanner. self.epi_series.append(info['series']) self.entry_map['epi'].append(entry) if not os.path.isdir(path): tmp_path = os.path.dirname(path) else: tmp_path = path self.epirt_paths.append(tmp_path) if self.fsl_flip: info['filetype'] = 'brik' else: info['filetype'] = self.tmplt['epi_file_format'] return OK def _GetImageInfo(self,path): hd = Header(path, scan=True) hdr = hd.hdr self.hdr = hdr if hdr is None: # Either a ref.dat file or it isn't an imaging file. if 'ref' in path and 'dat' in path: self.refdats[os.path.realpath(path)] = True info = {'type':'refdat'} return info else: return None elif hdr['filetype'] == 'dicom' and not path.endswith('.yaml'): # Write a yaml file to the raw data directory if possible. dirname, outfile = self._yaml_filename(path) yaml_name = '%s/%s' % (dirname, outfile) if not os.path.exists(yaml_name): # Create yaml file using dirname, # e.g., ../anatomicals/S2_EFGRE3D/s2_efgre3d.yaml try: hd.write_hdr_to_yaml('%s/%s' % (dirname,outfile)) except IOError: # This is a nonessential function, so ignore exceptions # such as access violations. pass elif hdr['filetype'] == 'dicom' or hdr['filetype'] == 'ge_ifile': if not os.path.isdir(path): path = os.path.dirname(path) shdr = hdr['subhdr'] nhdr = hdr['native_header'] self.shdr = shdr if 'dti' in shdr.get('PulseSequenceName','').lower() or hdr['mdim'] > 4 \ or 'dti' in nhdr.get('PulseSequenceFile',''): psdname = 'dti' else: psdname = os.path.basename((shdr.get('PulseSequenceName','').strip()).lower()) info = {'psdname':psdname, \ 'acqtime':shdr['AcqTime'], \ 'series':int(shdr['SeriesNumber']), \ 'plane':hdr['plane'].strip(), \ 'type':self.type.get(psdname,None), \ 'plane':hdr['plane'], \ 'acqtime':shdr['SeriesTime'], \ # 'fmapdir':None, \ 'refdat':None, \ 'imgfile':None, \ 'base':None, \ 'tdim':int(hdr['tdim']), \ 'echo_spacing':None, \ 'filetype':'brik', \ 'suffix':self.suffix.get(hdr['filetype'], 'brik'), \ 'data_filetype':hdr['filetype']} if info['type'] == 'localizer': # Don't process the localizer. return info if isinstance(info['acqtime'], int): info['acquisition_time'] = time.ctime(info['acqtime']) stat = apply( self.GetInfoMethods.get(info['type'], self._NoInfo), \ [info, path]) if stat: info = {'type':'break'} return info info['suffix'] = self.suffix.get(info['filetype'], 'brik') return info def StripSuffix(self, fname): if fname.endswith('.nii'): return fname[:-4] elif 'orig' in fname[-9:-5]: return fname[:-10] else: return fname def _GetFmapInfo(self): # Pair up each epi with a fieldmap. for epi in self.pfiles + self.epirt_paths: self.info[epi]['fmapname'] = None self.info[epi]['fmapdata'] = None self.info[epi]['anat_entry'] = self.anatref_entry if self.no_fmapcor: # Do not process fieldmaps if correction is to be skipped. continue for fmap in self.fmaps.keys(): if self.info[fmap]['plane'] == self.info[epi]['plane']: # Use the fieldmap acquired at the same plane. self.info[epi]['fmapname'] = '%s.nii' % self.fmaps[fmap] self.info[epi]['fmapdata'] = fmap self.info[epi]['anat_entry'] = fmap break else: for fmap in self.fmaps.keys(): # No fmap at same orientation, look for a sagittal if self.info[fmap]['plane'] == 'sagittal': self.info[epi]['fmapname'] = '%s.nii' % self.fmaps[fmap] self.info[epi]['fmapdata'] = fmap self.info[epi]['anat_entry'] = fmap break elif self.info[fmap]['plane'] == 'axial': # No sagittal unless this is overwritten. self.info[epi]['fmapname'] = '%s.nii' % self.fmaps[fmap] self.info[epi]['fmapdata'] = fmap self.info[epi]['anat_entry'] = fmap elif self.info[fmap]['plane'] == 'oblique': # No sagittal unless this is overwritten. self.info[epi]['fmapname'] = '%s.nii' % self.fmaps[fmap] self.info[epi]['fmapdata'] = fmap self.info[epi]['plane'] = 'oblique' self.info[epi]['anat_entry'] = fmap # Define which image the fieldmap should be normalized to. for entry in self.entry_map['fmap']: if self.info[entry]['anat_ref'] == 'self.anatomical' and \ self.anatomical is not None: self.info[entry]['anat_norm'] = self.StripSuffix(self.anatomical) else: self.info[entry]['anat_norm'] = self.info[entry]['anat_ref'] def _FindBaseEpi(self): # Find base image to be used for motion correction. for epi in self.pfiles + self.epirt_paths: anat_entry = None anat_ref = self.info[epi].get('anat_ref', None) if self.info[epi]['fmapdata'] is not None: # Found a fieldmap. Was it acquired before or after the data. anat_entry = self.info[epi]['fmapdata'] elif self.anat_ref: anat_entry = self.anat_ref else: for anat_entry in self.entry_map['anat']: if self.info[anat_entry].get('ref_name','') == anat_ref: # Found anat reference specified in template file. # This will overwrite the default value assigned in # _GetFmapInfo. break if anat_entry is None: raise RuntimeError('Could not find anatomical reference image.') self.info[epi]['anat_entry'] = anat_entry self.info[epi]['anat_matfile'] = self.info[anat_entry]['matfile'] epi_base = os.path.basename(self.info[epi]['imgfile_m']) self.info[epi]['matfile_m'] = '%s/%s.aff12.1D' % \ (self.info[epi]['outdir'], epi_base) self.info[epi]['matfile_mcat'] = '%s/%scat.aff12.1D' % \ (self.info[epi]['outdir'], epi_base) if self.info[anat_entry]['acqtime'] < self.info[epi]['acqtime'] \ or self.opts.base_firstepi: # Fieldmap acquired before EPI self.info[epi]['base'] = 'start' else: self.info[epi]['base'] = 'end' def _GetRefdat(self): for rfile in self.refdats.keys(): # Get times for ref.dat files with a time-stamp. words = rfile.replace('.','_').split('_') if len(words) == 6 and words[-2].count(':') == 20: # This file was time-stamped by the sequence. Get the # date and time. file name format: # ref_Sep_9_2007_11:28:32.dat rtime[rfile] = hms_to_secs(words[-2]) for pfile in self.pfiles: min_difftime = 1.e20 self.info[pfile]['refdat'] = None for rfile in self.refdats.keys(): if rfile[:3] == 'ref' and 'dat' in rfile: # This is a reference data file. First see if the orientation is # appended. If the file has neither a time-stamp nor a plane and # there is more than one ref.dat, the epi reconstruction will # be aborted. rinfo = {} ref_file = None if 'sag' in rfile and self.info[pfile]['plane'] == 'sagittal': # self.info[pfile]['refdat'] = rfile ref_file = rfile break elif 'cor' in rfile and self.info[pfile]['plane'] == 'coronal': # self.info[pfile]['refdat'] = rfile ref_file = rfile break elif 'axial' in rfile and self.info[pfile]['plane'] == 'axial': # self.info[pfile]['refdat'] = rfile ref_file = rfile break elif len(self.refdats.keys()) == 1: # Use the only one if that is all there is. ref_file = rfile epi_time = hms_to_secs(self.info[pfile]['acqtime'].split()[-2]) if epi_time - rtime[rfile] < min_difftime and \ rftime[rfile] > epi_time: # Use the reference file that acquired nearest to the EPI # but before it. min_difftime = epi_time - rtime[rfile] # self.info[pfile]['refdat'] = rfile ref_file = rfile if ref_file: # Found a candidate. if not self.info[pfile]['refdat']: # Haven't found one yet, use it. self.info[pfile]['refdat'] = ref_file else: # Found two. Choose one in the same directory. oldpath = os.path.dirname(self.info[pfile]['refdat']) newpath = os.path.dirname(ref_file) pfile_path = os.path.dirname(pfile) if oldpath == newpath: # Same path, use the old one. self.info[pfile]['refdat'] = ref_file elif newpath == pfile_path: self.info[pfile]['refdat'] = ref_file # else Do nothing, use existing choice. elif not os.path.exists(rfile): self.info[pfile]['refdat'] = None elif os.stat(rfile).st_size > 0: # This path is taken if no info is encoded in the file name. # Don't use empty ref.dat files. self.info[pfile]['refdat'] = rfile def _GetEpiOrder(self): # self.series_info = {} self.epi_series.sort() for series in self.epi_series: self.GetEpiAcqTimes(series) self.AssignEpiNames() def GetEpiAcqTimes(self, series): """ Fill structure for sorting acquisition times. """ # Find minimum and maximum start times for each acquistion in series. # min_time = None # max_time = None # sinfo = {} self.epi_times = {} for entry in self.entry_map['epi']: # Loop through each file in this series. if self.info[entry]['series'] == series and \ self.info[entry]['tdim'] > 2: # first_key = 'first_file_%s' % self.info[entry]['plane'] # last_key = 'last_file_%s' % self.info[entry]['plane'] # if min_time == None: # min_time = self.info[entry]['acqtime'] # max_time = self.info[entry]['acqtime'] # sinfo[first_key] = entry # sinfo[last_key] = entry # elif self.info[entry]['acqtime'] < min_time: # min_time = self.info[entry]['acqtime'] # sinfo[first_key] = entry # elif self.info[entry]['acqtime'] > min_time: # max_time = self.info[entry]['acqtime'] # sinfo[last_key] = entry # Relate each entry to its time of acquisition. self.epi_times[self.info[entry]['acqtime']] = entry # self.series_info[series] = sinfo # self.info[entry]['reconstruct'] = True # else: # self.info[entry]['reconstruct'] = False # return sinfo def AssignEpiNames(self): """ Assign names to each epi file based on information in the template. """ # Sort each run in the series by its acquisition time. epi_sort = self.epi_times.keys() epi_sort.sort() # Rewrite pfiles as an ordered list of p-files to be reconstructed. for idx in xrange(len(epi_sort)): entry = self.epi_times[epi_sort[idx]] info = self.info[entry] if info['data_filetype'] == 'ge_data': self.pfiles_recon.append(entry) info['run'] = '%0d' % (self.n_epi) self.n_epi = self.n_epi + 1 plane = info['plane'] if not self.epinames.has_key(plane): plane = 'any' n_epi = self.epinames[plane]['n_epi'] if n_epi > len(self.epinames[plane]['names'])-1: errstr = 'Not enough EPI names in template file' raise RuntimeError(errstr) filebase = os.path.basename(\ self.epinames[plane]['names'][n_epi]) epi_mf_outdir = os.path.dirname(\ self.epinames[plane]['names'][n_epi]) scratch_outdir = '/%s/%s' % \ (self.scratch, '/'.join(epi_mf_outdir.split('/')[2:])) tmp_outdir = '%s/%s' % \ (self.tmpdir, self.epinames[plane]['subdir'][n_epi]) epi_base = self.epinames[plane]['subdir'][n_epi] if self.keep_epi_raw: epi_r_outdir = scratch_outdir else: epi_r_outdir = tmp_outdir if self.keep_epi_mot: epi_m_outdir = scratch_outdir else: epi_m_outdir = tmp_outdir info['outdir'] = epi_mf_outdir if n_epi < len(self.epinames[plane]['names']): epiname = self.epinames[plane]['names'][n_epi] info['imgfile'] = '%s/%s' % (epi_r_outdir, filebase) else: info['imgfile'] = '%s/s%0d_epi_run%0d' % \ (epi_r_outdir, n_epi, idx+1) info['anat_ref'] = self.epinames[plane]['anat_ref'][n_epi] self.epinames[plane]['n_epi'] += 1 info['mot_file'] = '%s/%s_mtn.txt' % (epi_mf_outdir, filebase) info['censor_prefix'] = '%s/%s' % (epi_mf_outdir, filebase) info['imgfile_t'] = '%s/%s_t' % (epi_m_outdir, filebase) info['imgfile_m'] = '%s/%s_m' % (epi_m_outdir, filebase) info['imgfile_mf'] = '%s/%s_mf' % (epi_mf_outdir, filebase) info['skip'] = self.skip info['motion_ref_frame'] = self.tmplt['motion_ref_frame'] info['motion_interp'] = self.tmplt['epi_motion_interp'] if not info['motion_interp'].startswith('-'): info['motion_interp'] = '-%s' % info['motion_interp'] info['filetype'] = self.tmplt['epi_file_format'] self.info[entry] = info def DetermineBaseFrame(self): """ Determine base frame for motion correction. """ order = {\ 'axial':{'first_t':None, 'last_t':None, 'first_epi':None, 'last_epi':None},\ 'sagittal':{'first_t':None, 'last_t':None, 'first_epi':None, 'last_epi':None},\ 'coronal':{'first_t':None, 'last_t':None, 'first_epi':None, 'last_epi':None}} for entry in self.entry_map['epi']: plane = self.info[entry]['plane'] if order[plane]['first_t'] is None \ or self.info[entry]['acqtime'] < order[plane]['first_t']: order[plane]['first_t'] = self.info[entry]['acqtime'] order[plane]['first_epi'] = entry if order[plane]['last_t'] is None \ or self.info[entry]['acqtime'] > order[plane]['last_t']: order[plane]['last_t'] = self.info[entry]['acqtime'] order[plane]['last_epi'] = entry for entry in self.entry_map['epi']: # Determine base frame for motion correction. if self.info[entry]['type'] == 'epi' and \ self.info[entry]['tdim'] > 2: if self.opts.base_firstepi and len(self.entry_map['first_epi']) > 0: base_entry = self.entry_map['first_epi'][0] elif self.info[entry]['base'] == 'start': base_entry = order[self.info[entry]['plane']]['first_epi'] elif self.info[entry]['base'] == 'end': base_entry = order[self.info[entry]['plane']]['last_epi'] else: raise RuntimeError(\ 'Invalid value of base: __%s__' % self.info[entry]['base']) self.info[entry]['basefile'] = '%s' % \ (self.info[base_entry]['imgfile']) def DumpInfo(self): """ Dump the info dictionary to a yaml file. """ if self.logdir is None: return self.dumpfile = '%s/preprocess_info.yaml' % (self.logdir) try: f = open(self.dumpfile,'w') f.write(yaml.dump(self.info,default_flow_style=False, indent=4)) f.close() except IOError: self.errors = True errstr = 'Error accessing %s' % self.dumpfile raise IOError(errstr) self.LogErrors(errstr) def UnDumpInfo(self): """ Load the info dictionary from a yaml file. """ filename = '%s/preprocess_info.yaml' % self.logdir f = open(filename,'r') self.info = yaml.load(f.read()) f.close() class ProcData(DataInfo): def __init__(self, opts): self.dumpfile = None self.error_log = '' self.f_log = None self.f_crash = None self.f_bash = None self.opts = opts self.keep_epi_raw = opts.keep_epi_raw self.keep_epi_mot = opts.keep_epi_mot self.anatref_entry = None self.template_type = 'default' self.no_fmapcor = opts.no_fmapcor self.term_mesg = 'No problems detected' self.template_file = opts.template_file def Initialize(self, topdir, redo=False, verbose=False, dry_run=False, \ skip=None, scratch='/scratch', no_email=False, \ template_file=None): # Create default file object for early crashes. self.logfile = sys.stderr # Set default permission to 0775 os.umask(UMASK_FILE) # Initialize data structures. DataInfo.__init__(self, topdir, redo, verbose, skip, dry_run, scratch, \ no_email) if skip is None: self.skip = self.tmplt['skip'] self.redo = redo self.verbose = verbose self.dry_run = dry_run self.verbose = verbose self.no_email = no_email self.template_file = template_file # Create output directories. self.CreateDirs() self.DumpInfo() # Open the logfiles. self.logfile = '%s/preprocess.log' % self.logdir self.f_log = open(self.logfile,'w') self.f_log.write(self.verstring) self.f_log.write('# Written on %s\n\n' % \ datetime.today().strftime('%a %b %d, %Y; %X')) self.f_crash = open('%s/preprocess_failed.log' % self.logdir,'w') self.f_crash.seek(0,2) self.f_crash.write('Last written on %s\n\n' % \ datetime.today().strftime('%a %b %d, %Y; %X')) self.f_bash = open('%s/preprocess.bsh' % self.logdir,'w') self.f_bash.seek(0,2) self.f_bash.write('#!/bin/bash\n\n') self.f_bash.write('# Written on %s\n\n' % \ datetime.today().strftime('%a %b %d, %Y; %X')) def MakeDir(self, dirname): """ Create directory or exit on error. """ if os.path.exists(dirname): return try: os.umask(UMASK_DIR) os.makedirs(dirname) except OSError: self.errors = True errstr = '\nCould not create directory: %s ... ' % dirname self.LogErrors(errstr) raise OSError(errstr) os.umask(UMASK_FILE) def CreateDirs(self): # Create output directories if they don't already exist. # if not os.path.exists(self.procdir): # if self.verbose: # print "mkdir %s" % self.procdir # self.MakeDir(self.procdir) # Creat directory for initial epis # self.MakeDir(self.episetup_dir) # ltop = len(self.procdir) # dnames = [self.procdir, self.episetup_dir]] # First, create a list of directories. dnames = [] tags = ['', '_m', '_mf'] for entry in self.info.keys(): if self.info[entry]['type'] == 'epi': for tag in tags: fname = self.info[entry].get('imgfile%s' % tag, None) if fname is not None: dnames.append(os.path.dirname(fname)) else: if self.info[entry].get('outdir',None) is not None: dnames.append(self.info[entry]['outdir']) # Create them if they don't already exist. for dname in dnames: if not os.path.exists(dname): self.MakeDir(dname) if self.verbose: print 'mkdir %s' % dname def ExecCmd(self,cmd): self.f_bash.write("%s\n"%cmd) self.f_bash.flush() if not self.dry_run: execCmd(cmd, self.f_log, self.f_crash, self.verbose) def ConvertAnat(self): # Convert anatomical images to briks. if self.verbose: print 'Convert T1 and T2 images...' for entry in self.info: info = self.info[entry] if self.info[entry]['type'] in self.anat_types: key = self.info[entry]['type'] imgfile = self.info[entry]['imgfile'] cmd = 'convert_file %s %s %s %s' % (self.flip_opts, entry, \ imgfile, self.info[entry]['filetype']) checkfile = '%s%s' % (imgfile, self.info[entry]['suffix']) self.CheckExec(cmd, [checkfile]) if self.info[entry]['imgfile'] == \ self.info[entry]['anat_norm'] and \ self.skull_strip: cmd = "3dSkullStrip -input %s -prefix %s" % \ (checkfile, self.info[entry]['imgfile_skstrip']) checkfile = '%s+orig.BRIK' % \ (self.info[entry]['imgfile_skstrip']) self.CheckExec(cmd, [checkfile]) # def AlignAnat(self): for entry in self.entry_map['anat'] + self.entry_map['fmap']: if self.anatref_entry == None: # Nothing to align to. return if ((self.info[entry]['anat_norm'] != \ self.info[entry]['imgfile']) and \ (self.info[entry]['psdname'] == 'efgre3d') or \ self.info[entry]['psdname'] == '2dfast'): base_entry = self.anatref_entry basename = '%s%s' % (self.info[base_entry]['imgfile'], \ self.infile_suffix[self.info[base_entry]['filetype']]) # This is not the high res anatomical, but is either # fieldmap data or T1 IR data, so align to it. if self.info[entry]['psdname'] == '2dfast': fname = '%s' % self.info[entry]['magfile'] else: fname = self.info[entry]['imgfile'] infile = '%s%s' % \ (fname, self.infile_suffix[self.info[entry]['filetype']]) if os.path.exists(infile): cmd = "3dAllineate -prefix NULL -1Dmatrix_save %s " % \ self.info[entry]['matfile'] + \ "-base %s -source %s[0]" % (basename, infile) # 3dAllineate must create bricks if they are to store # transformation matrices. checkfile = '%s' % (self.info[entry]['matfile']) self.CheckExec(cmd, [checkfile]) else: # Assumed to be the same, write and identity transform. smat = "" for i in xrange(16): if i % 5: smat = '%s%14.9f' % (smat, 0.) else: smat = '%s%14.9f' % (smat, 1.) f = open('%s'%self.info[entry]['matfile'], 'w') f.write(smat) f.close() def ProcessDTI(self): # Convert anatomical images to briks. for entry in self.info: if self.info[entry]['type'] == 'dti': if self.verbose: print 'Processing DTI data in %s' % os.path.basename(entry) # dtiname = '%s/s%s_dti' % \ # (self.info[entry]['outdir'],self.info[entry]['series']) cmd = 'convert_file %s %s %s' % (entry, \ self.info[entry]['imgfile'], self.info[entry]['filetype']) fname = '%s%s' % \ (self.info[entry]['imgfile'], self.info[entry]['suffix']) self.CheckExec(cmd, [fname]) def MakeFieldmaps(self): # Create the fieldmap(s). if self.verbose: print 'Compute fieldmaps.' for entry in self.info: if self.info[entry]['type'] == 'fmap': if self.info[entry]['imgfile'] == None: # Fieldmap data not found. return # Make a magnitude image for use in checking registration. cmd = 'convert_file -f0 -m0 %s %s nii' % \ (entry, self.info[entry]['magfile']) self.CheckExec(cmd, [self.info[entry]['magfile'] + '.nii']) # Make fieldmap. Use separate loop in case make_fmap aborts. for entry in self.info: if self.info[entry]['type'] == 'fmap': fmapname = '%s' % self.info[entry]['imgfile'] if not os.path.exists('%s.nii' % fmapname) or self.redo: # Couldn't find or existing fmap, compute a new one. if self.verbose: extra_args = '-v' else: extra_args = '' cmd = 'make_fmap %s %s %s' % (extra_args, entry, fmapname) self.ExecCmd(cmd) def LinkAnat(self, srcdir): """Create link to structural image if it doesn't already exist.""" if self.anatomical is None: return self.LinkFiles(srcdir, self.anatomical) def LinkFiles(self, srcdir, target): if target.endswith('.BRIK'): linkfiles = [target, '%sHEAD' % target[:-4]] elif target.endswith('.HEAD'): linkfiles = [target, '%sHEAD' % target[:-4]] elif target.endswith('+orig'): linkfiles = ['%s.HEAD' % target, '%s.BRIK' % target] else: linkfiles = [target] for linkfile in linkfiles: linkname = '%s/%s' % (srcdir, os.path.basename(linkfile)) rel_linkdir = abspath_to_relpath(os.path.dirname(target), srcdir) rel_linkfile = '%s/%s' % (rel_linkdir, os.path.basename(linkfile)) if not os.path.exists(linkname) and not os.path.islink(linkname): cmd = 'cd %s && ln -s %s %s' % (srcdir, rel_linkfile, linkname) self.ExecCmd(cmd) def ExtractFirstEpi(self): # Extract the initial EPIs stored in dicom format. for entry in self.info: if self.info[entry]['type'] == 'first_epi': epiname = self.info[entry]['imgfile'] cmd = 'convert_file %s -f0 %s %s %s' % \ (self.flip_opts, entry,epiname, self.info[entry]['filetype']) fname = '%s%s' % (epiname, self.info[entry]['suffix']) self.CheckExec(cmd, [fname]) self.info[entry]['imgfile'] = fname # Create link to T1High self.LinkAnat(self.episetup_dir) def ReconEpis(self): # Reconstruct the EPIs. run = zeros(100) if self.verbose: print 'Reconstruct EPIs' for pfile in self.pfiles_recon: if self.info[pfile]['refdat'] is None: # Find the ref.dat file later. continue if self.info[pfile]['compression'] is not None: # Data are compressed, copy to tmp. compression = self.info[pfile]['compression'] pfile_decomp = '%s/%s' % (self.tmpdir, \ os.path.basename(self.info[pfile]['pfile_decomp'])) if os.path.exists(pfile_decomp): errstr = 'Attempting to overwrite existing p-file (%s)' % pfile_decomp + \ ' in ReconEpis' cmd = '%s %s > %s' % \ (decompress_cmds[compression], pfile, pfile_decomp) self.ExecCmd(cmd) else: # Create a link on /tmp to the pfile so the link to ref.dat will also # be on /tmp, (which is always writeable.) pfile_decomp = '%s/%s' % (self.tmpdir, os.path.basename(pfile)) if not os.path.exists(pfile_decomp): os.symlink(pfile, pfile_decomp) refname, refcmpress = self.CheckCompression( \ self.info[pfile]['refdat']) if refcmpress is not None: refdat_decomp = '%s/%s' % (self.tmpdir, os.path.basename(refname)) cmd = '%s %s > %s' % \ (decompress_cmds[refcmpress], \ self.info[pfile]['refdat'], refdat_decomp) self.ExecCmd(cmd) else: refdat_decomp = self.info[pfile]['refdat'] if refdat_decomp is not None: if refdat_decomp != 'ref.dat': # Create link bearing the file name epirecon_ex expects. refdat_link = '%s/ref.dat' % self.tmpdir if not os.path.exists(refdat_link): if self.verbose: print 'ln -s %s %s' % (refdat_decomp, refdat_link) if os.path.islink(refdat_link): # ref.dat is a broken symbolic link. if self.verbose: print 'rm %s' % ref_file os.remove(refdat_link) try: os.symlink(refdat_decomp, refdat_link) except OSError: self.errors = True pfile_link = '%s/%s' % (self.tmpdir, os.path.basename(pfile_decomp)) os.symlink(pfile_decomp, pfile_link) os.symlink(refdat_decomp, '%s/ref.dat' % self.tmpdir) series = int(self.info[pfile]['series']) run[series] = run[series] + 1 epiname = self.info[pfile]['imgfile'] cmd = 'epirecon_ex -F -f %s -NAME %s -fmt brik -skip %d' % \ (pfile_decomp,epiname,self.skip) fname = '%s+orig.BRIK' % epiname self.CheckExec(cmd, [fname]) else: errstr = '*******************************************\n' + \ 'No ref.dat file exists for %s\n' % pfile + \ '*******************************************\n' self.error_log = self.error_log + errstr self.f_crash.write(errstr) def ConvertRtEpis(self): # Convert epis reconstructed on the scanner. if self.verbose: print 'Convert EPIs to brik' for entry in self.entry_map['epi']: if (self.info[entry]['psdname'] == 'epirt' or \ self.info[entry]['psdname'] == 'epi') and \ self.info[entry]['data_filetype'] == 'dicom': series = self.info[entry]['series'] if self.info[entry]['skip'] > 0: skip = '--skip=%s' % self.info[entry]['skip'] else: skip = '' cmd = 'convert_file %s %s %s brik' % \ (skip, \ entry, self.info[entry]['imgfile']) # self.info[entry]['filetype']) checkname = '%s+orig.BRIK' % (self.info[entry]['imgfile']) # self.suffix[self.info[entry]['filetype']]) self.CheckExec(cmd, [checkname]) def CorrectMotion(self): # Motion correction if self.verbose: print "Correct for motion" # for entry in self.info: for entry in self.entry_map['epi']: if self.info[entry]['type'] == 'epi' and \ self.info[entry].has_key('imgfile_m'): # Always use brik for 3dDeconvolve. ending = '+orig' infile = '%s%s' % (self.info[entry]['imgfile'], ending) prefix = self.info[entry]['imgfile_m'] if self.info[entry]['base'] == 'start': base = self.info[entry]['motion_ref_frame'] else: base = self.info[entry]['tdim'] - self.info[entry]['skip']-1 base = ('%d' % base).replace(' ','') anat_entry = self.info[entry]['anat_entry'] if anat_entry is not None and \ (self.info[anat_entry].get('anat_norm', '') != \ self.info[anat_entry].get('imgfile', '')) \ and self.epi_align: # Include additonal transformation in motion correction. # First, time-shift and compute matrices. cmd = '3dvolreg -prefix NULL -1Dmatrix_save %s -twopass ' %\ ( self.info[entry]['matfile_m']) + \ '-tshift -verbose -base %s%s[%s] -dfile %s %s' % \ (self.info[entry]['basefile'], ending, base, \ self.info[entry]['mot_file'], infile) # self.infile_suffix[self.info[entry]['filetype']], \ # Compute the transformation matrices. self.CheckExec(cmd, [self.info[entry]['matfile_m']]) # Catenate with transformation from base to anatomical. cmd = 'cat_matvec -ONELINE %s -I %s -I > %s' %\ (self.info[anat_entry]['matfile'], \ self.info[entry]['matfile_m'], \ self.info[entry]['matfile_mcat']) self.CheckExec(cmd, [self.info[entry]['matfile_mcat']]) # Time-shift the EPIs if self.skip == 0: ignore = 3 else: ignore = self.skip cmd = '3dTshift -prefix %s -ignore %d %s' % \ (self.info[entry]['imgfile_t'], ignore, infile) checknames = ['%s+orig.BRIK'%self.info[entry]['imgfile_t'],\ '%s+orig.HEAD'%self.info[entry]['imgfile_t']] self.CheckExec(cmd, checknames) # Realign. cmd = '3dAllineate -prefix %s -interp cubic ' % prefix + \ '-1Dmatrix_apply %s -base %s%s[%s] %s+orig' % \ (self.info[entry]['matfile_mcat'], \ self.info[entry]['basefile'], ending, base, \ self.info[entry]['imgfile_t']) self.CheckExec(cmd, ['%s+orig.BRIK'%prefix, \ '%s+orig.HEAD'%prefix]) else: cmd = '3dvolreg -prefix %s -twopass -tshift %s -verbose ' %\ (prefix, self.info[entry]['motion_interp']) + \ '-base %s%s[%s] -dfile %s %s' % \ (self.info[entry]['basefile'], ending, \ base, self.info[entry]['mot_file'], infile) fname = '%s+orig.BRIK' % prefix self.CheckExec(cmd, ['%s+orig.BRIK' % prefix, \ '%s+orig.HEAD' % prefix]) # Create link to T1High from each unique directory name. dnames = {} dnames[os.path.dirname(self.info[entry]['imgfile'])] = True dnames[os.path.dirname(self.info[entry]['imgfile_m'])] = True dnames[os.path.dirname(self.info[entry]['imgfile_mf'])] = True for dname in dnames.keys(): self.LinkAnat(dname) if not self.info[entry].get('fmapname', None): if self.fsl_flip: # Must be nifti prefix = '%s/%s' % \ (os.path.dirname(self.info[entry]['imgfile_mf']), \ os.path.basename(self.info[entry]['imgfile_m'])) self.FSLFlip(self.info[entry]['imgfile_m'], prefix) elif not self.keep_epi_mot: # Copy motion-corrected images from /tmp to output directory infile = '%s+orig.BRIK' % (self.info[entry]['imgfile_m']) outfile = '%s/%s+orig.BRIK' % ( \ os.path.dirname(self.info[entry]['imgfile_mf']),\ os.path.basename(self.info[entry]['imgfile_m'])) cmd = "cp %s %s" % (infile, outfile) self.CheckExec(cmd, [outfile]) infile = '%s+orig.HEAD' % (self.info[entry]['imgfile_m']) outfile = '%s/%s+orig.HEAD' % ( \ os.path.dirname(self.info[entry]['imgfile_mf']),\ os.path.basename(self.info[entry]['imgfile_m'])) cmd = "cp %s %s" % (infile, outfile) self.CheckExec(cmd, [outfile]) # Create censor file for use in 3dDeconvolve. self.JumpCensor(entry) def JumpCensor(self, entry): if self.censor_interleave: interleave = '--interleave' input_file = '%s+orig' % self.info[entry]['imgfile'] else: interleave = 'interleave' input_file = self.info[entry]['mot_file'] cmd = "jump_censor --prefix=%s %s --store-plot --threshold=%f %s" % \ (self.info[entry]['censor_prefix'], \ interleave, \ self.censor_thresh, \ input_file) self.CheckExec(cmd, ['%s_censor.1D' % self.info[entry]['censor_prefix']]) def CheckExec(self, cmd, checkname): """ Check if output file exists, then execute commmand. If there is more than one output file, the command will be executed if at least one is missing. """ gone = False names = [] for name in checkname: if '+orig' in name: if name.endswith('+orig'): names.append('%s.HEAD' % name) names.append('%s.BRIK' % name) elif name.endswith('HEAD'): names.append(name) newname = name[:-4] + 'BRIK' if newname not in checkname: names.append(newname) elif name.endswith('BRIK'): newname = name[:-4] + 'HEAD' if newname not in checkname: names.append(newname) names.append(name) else: names.append(name) for name in names: if not os.path.exists(name): gone = True elif self.redo: os.remove(name) if self.redo or gone: self.ExecCmd(cmd) if '+orig.' in checkname[0]: append_history_note(checkname[0][:-5], cmd) def FieldmapCorrection(self): if self.verbose: print "Fieldmap correction" for entry in self.entry_map['epi']: if self.info[entry].get('fmapname', None) is None: # No fieldmap, save motion corrected image. infile = self.info[entry]['imgfile_m'] prefix = '%s/%s' % \ (self.info[entry]['outdir'], \ os.path.basename(self.info[entry]['imgfile_m'])) if self.fsl_flip: self.FSLFlip(infile, prefix) else: outfile = '%s+orig.BRIK' % prefix cmd = 'cp %s+orig.BRIK %s' % (infile, outfile) self.CheckExec(cmd, outfile) outfile = '%s+orig.HEAD' % prefix cmd = 'cp %s+orig.HEAD %s+orig.HEAD' % (infile, prefix) self.CheckExec(cmd, outfile) elif self.info[entry]['type'] == 'epi': infile = '%s+orig' % (self.info[entry]['imgfile_m']) fmapdir = self.info[entry]['outdir'] prefix = self.info[entry]['imgfile_mf'] cmd = 'fieldmap_correction -m --tag=f -t%s %s %s %s %s' % \ (self.info[entry]['filetype'], \ self.info[entry]['fmapname'], \ self.info[entry]['echo_spacing'], \ os.path.dirname(prefix), \ infile) if self.info[entry]['filetype'] == 'brik': fname = '%s%s' % (prefix, self.info[entry]['suffix']) else: fname = '%s.nii' % prefix if self.redo: self.ExecCmd(cmd) elif not os.path.exists(fname): self.ExecCmd(cmd) if self.fsl_flip: self.FSLFlip(prefix, prefix) # Create links from outut directory to files on /scratch for tgt in (self.info[entry]['imgfile'],self.info[entry]['imgfile_m']): if '/scratch' in tgt: self.LinkFiles(self.info[entry]['outdir'], \ tgt+self.info[entry]['suffix']) def FSLFlip(self, infile, prefix): # Flip axes to orientation fslview expects. cmd = '3dresample -orient LPI -prefix %s.nii -inset %s+orig' % \ (prefix, infile) self.CheckExec(cmd, ['%s.nii' % prefix]) fname = '%s+orig.BRIK' % infile if os.path.exists(fname): os.remove(fname) fname = '%s+orig.HEAD' % infile if os.path.exists(fname): os.remove(fname) def Chown(self): # Change ownership to group read read/write. cmd = 'chmod -R 0775 %s' % self.procdir self.ExecCmd(cmd) def LogErrors(self, errstr): self.error_log = self.error_log + errstr if self.opts.verbose: sys.stderr.write('%s\n' % errstr) if self.f_crash is not None: self.f_crash.write('\n%s\n' % errstr) if self.f_log is not None: self.f_log.write('\n%s\n' % errstr) def LogProcess(self): time = datetime.today().strftime('%a %Y%b%d %X') # Get user name. f = os.popen("whoami","r") user = f.read().strip() f.close() # Read version from the release directory. # f = open(c.VERSION_FILE,'r') # self.version = f.read().strip() # f.close() entry = '%s\t%s\t%s\t%s\n' % (time, self.topdir, user, self.version) # Append info to the exams file. try: f = open(c.EXAMS_FILE,'a+') f.seek(0, 2) f.write(entry) f.close() except: pass def CleanUp(self): if (not self.keep_epi_raw or not self.keep_epi_mot) \ and not self.opts.debug_tmp: cmd = "rm -rf %s" % self.tmpdir if self.opts.verbose: print cmd os.system(cmd) overall_msg = self.SummaryErrorMessage() if self.tmplt and not self.no_email: EmailResults(self.tmplt['email'], overall_msg, \ self.topdir, self.dumpfile, self.logfile) # Write the error message to the log file. if self.f_log is None: # Log file not opened yet, do it now. if self.logdir is not None: logfile = '%s/preprocess.log' % self.logdir f_log = open(logfile,'w') f_log.write('\n%s\n' % overall_msg) f_log.close() else: self.f_log.write('\n%s\n' % overall_msg) sys.exit() def SummaryErrorMessage(self, error_log=None): if error_log is None: error_log = self.error_log server = socket.gethostname().split('.')[0] mssg = '\nPreprocessing script complete for data in %s\n\nServer: %s\n'\ % (self.topdir, server) # Log time. ms = time.time() ms = int(1000*(ms - int(ms))) mssg += '\nTime: %s:%03d\n' % \ (datetime.today().strftime('%a %b %d, %Y; %X'), ms) if len(error_log) > 0: mssg += 'Command: %s\n\nSummary:\n' % (' '.join(sys.argv)) lines = error_log.split('\n') for line in lines: if line.startswith('Description:'): mssg += line[12:] mssg += '\n\nDetails:' + error_log else: mssg += '\nNo problems detected (this does NOT imply that everything was computed.).\n\n' return mssg def key_callback(option,opt_str,value,parser): """ Used with optparser for multiple arguments of the same type. """ if "--epi-key" in opt_str: parser.values.epi_keys.append(value) elif "--exclude" in opt_str: parser.values.exclude_paths.append(value) def main(user_email): proc = None usage = "preprocess path Enter --help for a list of options." + \ "Default action is to process everything." optparser = OptionParser(usage) optparser.add_option( "-v", "--verbose", action="store_true", \ dest="verbose",default=False, \ help='Print stuff to screen.') optparser.add_option( "", "--no-email", action="store_true", \ dest="no_email",default=False, \ help="Don't send any emails.") optparser.add_option( "-R", "--recompute", action="store_true", \ dest="redo",default=False, help='Recompute all files.') optparser.add_option("", "--exclude", action="callback", \ callback=key_callback, dest="exclude_paths", type="string", \ default=[], help=\ 'A path that is to be excluded from the analysis. Both ' + \ 'relative and absolute paths can be entered. If EXCLUDE_PATH is ' +\ 'a directory, all files in the directory and in its subdirectories '+\ 'will be excluded.') optparser.add_option("", "--epi-key", action="callback", \ callback=key_callback, dest="epi_keys", type="string", \ default=[], help=\ 'Keyword and value that must be present in the p-file header if it' + \ 'is to be processed. For example, suppose epis were collected with '+\ 'for two paradigms during the same session. If one contained 198 ' + \ 'frames and the other 210, specfying the keyword tdim:198 would '+\ 'cause the program to only process EPI runs containing 198 frames' + \ 'Valid keywords are "tdim" (number of frames), "SeriesNumber", '+\ 'and/or "plane" (axial, sagittal, or coronal). Multiple keywords' + \ 'can be entered (e.g., "--epi_key=tdim:198 --epi-key=plane:sagittal")') optparser.add_option( "", "--hog-disk", action="store_true", \ dest="hog",default=False, \ help='Save all intermediate files. Same as --keep-epi-all.') optparser.add_option( "", "--keep-epi-raw", action="store_true", \ dest="keep_epi_raw",default=False, \ help='Keep uncorrected EPI images. They will be stored ' + \ 'in a shadow directory on /scratch, e.g., if output ' + \ 'directory is /study/mystudy/processed/subj001, images ' + \ 'will be stored in /scratch/mystudy/processed/subj001') optparser.add_option( "", "--keep-epi-mot", action="store_true", \ dest="keep_epi_mot",default=False, \ help='Keep motion-corrected (but not yet fieldmap-' +\ 'corrected EPI images. They will be stored ' + \ 'in a shadow directory on /scratch, e.g., if output ' + \ 'directory is /study/mystudy/processed/subj001, images ' + \ 'will be stored in /scratch/mystudy/processed/subj001') optparser.add_option( "", "--keep-epi-all", action="store_true", \ dest="keep_epi_all",default=False, \ help='Keep all EPI images, including intermediate images ' + \ 'Intermediate EPIs will be stored ' + \ 'in a shadow directory on /scratch, e.g., if output ' + \ 'directory is /study/mystudy/processed/subj001, images ' + \ 'will be stored in /scratch/mystudy/processed/subj001') # optparser.add_option( "-o", "--out_dir", action="store", dest="outdir", \ # type="string",default=None,help="Directory where output "+ \ # "files are to be written. Defaults to current working directory.") optparser.add_option( "", "--base-firstepi", action="store_true", \ dest="base_firstepi", default=False, help=\ 'Align EPI runs to reference EPI image instead of first or last '+\ 'timeseries image.' ) optparser.add_option( "", "--dry_run", action="store_true", \ dest="dry_run",default=False, help= \ "Analyze data and create output yaml file but done't compute anything.") optparser.add_option( "", "--debug-tmp", action="store_true", \ dest="debug_tmp",default=None, help=\ "Write tmp files to /tmp/debug_tmp. For debugging purposes only.") optparser.add_option( "", "--skull-strip", action="store_true", \ dest="skull_strip",default=False, \ help='Skull strip the anatomical reference image.') optparser.add_option( "", "--align-epis", action="store_true", \ dest="align_epis",default=False, \ help='Align EPIs with anatomical reference by concatenating' + \ 'transformation computed for motion correction with the ' + \ 'transformation from the epis "anat_ref" image to the ' + \ 'high-res anatomical.') optparser.add_option( "-s", "--skip", action="store", dest="skip", \ type="int",default=None,\ help="Number of frames to skip at the beginning of each run.") optparser.add_option( "", "--template", action="store", \ dest="template_file", type="str",default=None,\ help="The one and only template file to be used.") optparser.add_option( "", "--scratch", action="store", dest="scratch", \ type="str",default='/scratch',\ help="Name of scratch partition. Destination for " + \ "intermediate time-series data.") optparser.add_option( "-o", "--output_dir", action="store", dest="outdir", \ type="string",default=None,\ help="Directory where output data are written. " + \ "This overides the value in the template file.") optparser.add_option( "", "--followon-script", action="store", \ dest="followon_script", type="string",default=None,\ help="Script to be run after preprocess complets." + \ "It will be called with the output path as the only argument.") optparser.add_option( "", "--no_fmapcor", action="store_true", \ dest="no_fmapcor",default=False, \ help='Skip fieldmpap correction even if a fieldmap exists.') optparser.add_option( "-a", "--anat", action="store_true", \ dest="anat",default=False, \ help='Convert structural images.') optparser.add_option( "-d", "--dti", action="store_true", \ dest="dti",default=False, \ help='Process dti images.') optparser.add_option( "-e", "--epi", action="store_true", \ dest="epi",default=False, \ help='Reconstruct epi images.') optparser.add_option( "-m", "--motion", action="store_true", \ dest="motion",default=False, \ help='Motion-correct epi images.') optparser.add_option( "-f", "--fmap", action="store_true", \ dest="fmap",default=False, \ help='Fieldmap-correct epi images.') optparser.add_option( "-c", "--compute_fmap", action="store_true", \ dest="cfmap",default=False, \ help='Compute fieldmap-correction.') optparser.add_option( "-A", "--Anat", action="store_true", \ dest="Anat",default=False, \ help='Convert structural images, (recompute).') optparser.add_option( "-D", "--Dti", action="store_true", \ dest="Dti",default=False, \ help='Process dti images, (recompute).') optparser.add_option( "-E", "--Epi", action="store_true", \ dest="Epi",default=False, \ help='Reconstruct epi images, (recompute).') optparser.add_option( "-M", "--Motion", action="store_true", \ dest="Motion",default=False, \ help='Motion-correct epi images, (recompute).') optparser.add_option( "-F", "--Fmap", action="store_true", \ dest="Fmap",default=False, \ help='Fieldmap-correct epi images, (recompute).') optparser.add_option( "-C", "--Compute_fmap", action="store_true", \ dest="Cfmap",default=False, \ help='Compute fieldmap-correction, (recompute).') optparser.add_option( "-V", "--version", action="store_true", \ dest="show_version",default=None, help="Display svn version.") error_log = "" opts, args = optparser.parse_args() if opts.show_version: sys.stdout.write('%s\n' % ID) sys.exit() if len(args) != 1: # optparser.error( \ # "Expecting 1 argument:\nEnter 'preprocess --help' for usage." ) f = open(c.FMRI_TOOLS_VERSION_FILE,'r') version = f.read().strip() f.close() print "\npreprocess version %s" % version errstr = "Expecting 1 argument:\nEnter 'preprocess --help' for usage." error_log = error_log + errstr # proc.LogErrors(errstr) raise RuntimeError(errstr) # Name the directories. topdir = os.path.abspath(args[0]) print 'Processing %s' % topdir if not os.path.exists(topdir): errstr = 'preprocess: Cannot access data: %s\n' % topdir error_log = error_log + errstr sys.stderr.write(errstr) # proc.LogErrors(errstr) raise IOError(errstr) # Store the path so we can tell what was done today. f = open(c.exams_file,'r') # f.seek(0,2) lines = f.readlines() f.close() outpaths = [] outlines = [] for line in lines: words = line.strip().split() if words[-1] not in outpaths: outpaths.append(words[-1]) outlines.append(line) if topdir not in outpaths: timestamp = datetime.today().strftime('%a %b %d, %Y; %X') outlines.append('%s; %s\n' % (timestamp,topdir)) f = open(c.exams_file,'w') f.writelines(outlines) f.close() # Set flag to compute the item. opts.anat = opts.anat or opts.Anat opts.dti = opts.dti or opts.Dti opts.cfmap = opts.cfmap or opts.Cfmap opts.epi = opts.epi or opts.Epi opts.motion = opts.motion or opts.Motion opts.fmap = opts.fmap or opts.Fmap if opts.keep_epi_all or opts.hog: opts.keep_epi_raw = True opts.keep_epi_mot = True # Don't process dti data. opts.dti = False if not (opts.anat or opts.dti or opts.fmap or opts.epi or opts.motion or \ opts.cfmap): # If no options specified, do everything that hasn't been done.. opts.anat = True opts.dti = True opts.cfmap = True opts.epi = True opts.motion = True opts.fmap = True if opts.redo: # Recompute everything opts.Anat = opts.Dti = opts.Fmap = opts.Epi = \ opts.Motion = opts.Cfmap = True try: # Create processing object. proc = ProcData(opts) proc.error_log = error_log except RuntimeError, errmsg: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = 'preprocess: Error during initialization.\n%s\n' % errmsg sys.stderr.write(errstr) error_log = error_log + errstr sys.exit(1) try: # Initialize data strucures by scanning all subdirectories. proc.Initialize(topdir, opts.redo, opts.verbose, \ opts.dry_run, opts.skip, opts.scratch, opts.no_email, \ opts.template_file) if proc.template_type == 'default': sys.stderr.write('No template found. Exiting\n') proc.term_mesg = '\nAbnormal termination.\n' sys.exit(1) if opts.dry_run: # Test run. info structure was computed, so exit. proc.DumpInfo() sys.stdout.write(\ '\nFinished initialization. Information used to guide ' + \ 'preprocessing was\nwritten to %s\n' % proc.dumpfile) sys.exit() # user_email = proc.server_userid if opts.anat: proc.redo = opts.Anat | opts.redo proc.ConvertAnat() proc.redo = opts.redo if opts.dti: proc.redo = opts.Dti | opts.redo proc.ProcessDTI() proc.redo = opts.redo if opts.cfmap: proc.redo = opts.Cfmap | opts.redo proc.MakeFieldmaps() proc.redo = opts.redo if opts.epi: proc.redo = opts.Epi | opts.redo proc.ExtractFirstEpi() proc.AlignAnat() proc.ReconEpis() proc.ConvertRtEpis() proc.redo = opts.redo if opts.motion: proc.redo = opts.Motion | opts.redo proc.CorrectMotion() proc.redo = opts.redo if opts.fmap: proc.redo = opts.Fmap | opts.redo proc.FieldmapCorrection() proc.redo = opts.redo proc.CleanUp() except ExecError, err: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = except_msg(err.errmsg) proc.LogErrors(errstr) proc.CleanUp() sys.exit() except RuntimeError, errmsg: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = except_msg(errmsg) proc.LogErrors(errstr) proc.CleanUp() sys.exit() except ScannerError, err: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = except_msg('--- Error in template file. ---') proc.LogErrors(errstr) proc.CleanUp() # Fatal error, exit sys.exit(1) except SystemExit: pass # proc.term_mesg = '\nAbnormal termination.\n' except KeyboardInterrupt: proc.term_mesg = '\nTerminated by ctrl-c\n' sys.exit() except: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = except_msg(__name__) proc.LogErrors(errstr) proc.CleanUp() # Change group. try: if sys.platform != 'darwin' and hasattr(proc, 'procdir'): # Won't work on the apples. cmd = 'chgrp -R dusers %s' % proc.procdir execCmd(cmd,proc.f_log,proc.f_crash,proc.verbose) except (ExecError, RuntimeError), err: # This is a nonfatal error, log it and continue. proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' proc.LogErrors("Could not change group ID of output files to dusers") # Change permissions. try: uid = os.getuid() if hasattr(proc, 'procdir'): for root,dummy,files in os.walk(proc.procdir): for fname in files: fullname = '%s/%s' % (root, fname) if os.path.exists(fullname): st = os.stat(fullname) if st.st_uid == uid: # File ownership OK, add user and group write access. os.chmod(fullname, S_IRWXU | S_IRWXG ) except ExecError, err: # This is a nonfatal error. proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' proc.LogErrors("Could not change ownership of output files.") # Dump status information in the info dictionary to disk and compress # the raw data. try: proc.DumpInfo() except IOError, errstr: proc.errors = True proc.term_mesg = '\nAbnormal termination.\n' errstr = except_msg(errstr) proc.LogErrors(errstr) proc.CleanUp() # Tell them its done. if proc.tmplt is not None: overall_msg = proc.SummaryErrorMessage() if not proc.no_email: EmailResults(proc.tmplt['email'], overall_msg, \ proc.topdir, proc.dumpfile, proc.logfile) if proc.errors and not proc.no_email: # Mail me the errors. EmailResults('ollinger@wisc.edu', overall_msg, \ proc.topdir, proc.dumpfile, proc.logfile) else: if proc is None: print 'Abnormal termination.' else: print proc.term_mesg if opts.followon_script is not None: cmd = '%s %s' % (opts.followon_script, proc.procdir) print cmd os.system(cmd) if __name__ == '__main__': if 'linux' in sys.platform: libc = dl.open('/lib/libc.so.6') libc.call('prctl',15,'preprocess',0,0,0) try: user_email = None main(user_email) except RuntimeError, errmsg: sys.stderr.write('\nError occured in preprocess.\n%s\n\n' % errmsg) EmailResults(user_email, errmsg, sys.argv, None, '/dev/null') sys.exit(1)