#!/usr/bin/env python # system libraries import sys, os, glob # AFNI libraries import afni_util as UTIL import option_list as OL # ---------------------------------------------------------------------- # globals g_help_string = """ ============================================================================= compute_OC_weights.py - compute weights for optimally combining echoes Given echo times (in a text file) and one run of multi-echo EPI data, compute a dataset that can be used to combine the echoes. The weight dataset would have one volume per echo, which can be used to combine the echoes into a single dataset. The same echoes can be applied to all runs. 3dcalc -a echo1+tlrc -b echo2+tlrc -c echo3+tlrc \ -d weights+tlrc'[0]' -e weights+tlrc'[1]' -f weights+tlrc'[2]' \ -expr 'a*d+b*e+c*f' -prefix opt.combined ---------------------------------------------------------------------- These computations are based on the system of equations from the summer 2017 presentation by Javier Gonzalez-Castillo. After solving: log(mean(S(TE_1))) ~= -mean(R2s(x))*TE_1 + log(So(x)) log(mean(S(TE_2))) ~= -mean(R2s(x))*TE_2 + log(So(x)) log(mean(S(TE_3))) ~= -mean(R2s(x))*TE_3 + log(So(x)) then T2* = 1/mean(R2s(x)), and weights come from: TE_n*e^-(TE_n/T2*) w(TE_n) = ------------------------- sum_n[TE_n*e^-(TE_n/T2*)] Bad, naughty voxels are defined as those with either negative T2* values, or for which the sum of the weights is not sufficiently close to 1, which would probably mean that there were computational truncation errors, likely due to R2s being very close to 0. so "fail" if mean(R2s) <= 0 or abs(1-sum[weights]) > 'tolerance' In such cases, the weights will default to 1/number_of_echoes. ---------------------------------------------------------------------- ------------------------------------------ examples: ------------------------------------------ terminal options: -help : show this help -hist : show the revision history -ver : show the version number ------------------------------------------ required parameters: -echo_times_file FILE - specify file with echo times (e.g. it could contain 15 30.5 41) -echo_times_file echo_times.1D -echo_dsets D1 D2 D3 - specify one run of multi-echo EPI data, e.g.: -echo_dsets pb03.SUBJ.r01.e*.volreg+tlrc.HEAD ------------------------------------------ other options: -def_to_equal yes/no - specify whether to default to equal weights (default = yes) -def_to_equal no In the case where T2* seems huge or <= 0, or if the sum of the fractional weights is not close to 1 (see -tolerance), one might want to apply default weights equal to 1/num_echoes (so echoes are weighted equally). -prefix PREFIX - specify prefix of resulting OC weights dataset -prefix OC.weights.SUBJ Specify the prefix for the main output (weights) dataset. -t2_star_limit LIMIT - specify limit for T2* values (default 300) -t2_start_limit 100 This is probably not important. -work_dir RDIR - specify directory to compute results in -work_dir temp.OC -verb - increase verbosity of output --------------------------------------------------------------------------- R Reynolds, Feb. 2018 Thanks to Javier Gonzalez-Castillo ============================================================================= """ g_todo = """ todo list: """ g_history = """ compute_OC_weights.py history: 0.0 Feb 12, 2018 - started (tcsh script) 0.1 Feb 16, 2018 - initial version (translate tcsh to python) """ g_version = "compute_OC_weights.py version 0.1, Feb 16, 2018" class MyInterface: """interface class for MyLibrary (whatever that is) """ def __init__(self, verb=1): # main variables self.valid_opts = None self.user_opts = None # control self.verb = 1 # infile name parsing self.infiles = [] # initialize valid_opts self.valid_opts = self.get_valid_opts() def get_valid_opts(self): vopts = OL.OptionList('valid opts') # short, terminal arguments vopts.add_opt('-help', 0, [], helpstr='display program help') vopts.add_opt('-hist', 0, [], helpstr='display the modification history') vopts.add_opt('-ver', 0, [], helpstr='display the current version number') # require options vopts.add_opt('-echo_dsets', -1, [], helpstr='required: one run of multi-echo data') vopts.add_opt('-echo_times', -1, [], helpstr='required: list of echo times') vopts.add_opt('-prefix', 1, [], helpstr='prefix for main output: echo weight dataset') vopts.add_opt('-sum_weight_tolerance', 1, [], helpstr='tolerance for summed weights being close to 1.0') vopts.add_opt('-work_dir', 1, [], helpstr='directory to do work in') vopts.add_opt('-verb', 1, [], helpstr='set the verbose level (def=1)') vopts.sort() return vopts def process_options(self): """return 1 on valid and exit return 0 on valid and continue return -1 on invalid """ argv = sys.argv # process any optlist_ options self.valid_opts.check_special_opts(argv) # process terminal options without the option_list interface # (so that errors are not reported) # if no arguments are given, do default processing if '-help' in argv or len(argv) < 2: print g_help_string return 1 if '-hist' in argv: print g_history return 1 if '-show_valid_opts' in argv: self.valid_opts.show('', 1) return 1 if '-ver' in argv: print g_version return 1 # ============================================================ # read options specified by the user self.user_opts = OL.read_options(argv, self.valid_opts) uopts = self.user_opts # convenience variable if not uopts: return -1 # error condition # ------------------------------------------------------------ # process verb first val, err = uopts.get_type_opt(int, '-verb') if val != None and not err: self.verb = val # ------------------------------------------------------------ # process options sequentially, to make them like a script errs = 0 for opt in self.user_opts.olist: # check for anything to skip if opt.name == '-verb': pass elif opt.name == '-infiles': self.infiles, err = uopts.get_string_list('', opt=opt) if self.infiles == None or err: print '** failed to read -infiles list' errs +=1 # allow early and late error returns if errs: return -1 # ------------------------------------------------------------ # apply any trailing logic if len(self.infiles) < 1: print '** missing -infiles option' errs += 1 if errs: return -1 return 0 def do_stuff(self): """main function to process input """ print '-- have %d files to process' % len(self.infiles) errs = 0 # check file existence first for ifile in self.infiles: if ifile in ['-', 'stdin']: pass elif not os.path.isfile(ifile): print '** input file not found: %s' % ifile errs += 1 if errs: return 1 return 0 def main(): me = MyInterface() if not me: return 1 rv = me.process_options() if rv > 0: return 0 # exit with success if rv < 0: # exit with error status print '** failed to process options...' return 1 if me.do_stuff(): return 1 return 0 if __name__ == '__main__': sys.exit(main())