#!/usr/bin/env python
# -*- coding: latin-1 -*-
#
# long_mris_slopes
#
# script to fit within-subject slopes into longitudinal data
#
# Original Author: Martin Reuter
# CVS Revision Info:
# $Author: mreuter $
# $Date: 2012/05/30 22:57:25 $
# $Revision: 1.38 $
#
# Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
#
# Terms and conditions for use, reproduction, distribution and contribution
# are found in the 'FreeSurfer Software License Agreement' contained
# in the file 'LICENSE' found in the FreeSurfer distribution, and here:
#
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
#
# Reporting: freesurfer@nmr.mgh.harvard.edu
#
#
import warnings
warnings.filterwarnings('ignore', '.*negative int.*')
import os
import sys
import shlex
import optparse
import logging
import subprocess
import tempfile
import shutil
from subject_info import *
from LongQdecTable import *
# logging
ch = logging.StreamHandler()
#create logger
slopelogger = logging.getLogger("long_mris_slopes")
slopelogger.setLevel(logging.INFO)
slopelogger.addHandler(ch)
ARGTEXT = """
REQUIRED ARGUMENTS
--qdec <name> qdec.table.dat file with first columns: fsid fsid-base
--meas <name> Input curv file, e.g. thickness (except if only --do-label)
--hemi ?h Hemisphere (lh or rh)
One or more of the following:
--do-avg Compute and output the temporal average (recommended)
--do-rate Compute and output the rate (recommended)
--do-pc1fit Compute and output the pct. change (w.r.t. tp1 computed from line fit)
--do-pc1 Compute and output the pct. change (w.r.t. tp1)
--do-spc Compute and output the sym. pct. change (w.r.t. temp. average) (recommended)
--do-label Compute and output intersected cortex label (recommended)
OPTIONAL ARGUMENTS
--do-stack Save the stacked within subject file (time series) (recommended)
--fwhm <int> Smooth the data (recommended for percent change maps)
--nosmooth Do not smooth (to overwrite smoothing requirement for pc1 or spc)
--time <name> Variable name of time variable (e.g. age)
--in-label <name> Use pre-existing label for smoothing (default: intersect cortex labels)
--qcache <qtarget> Create cache for qdec (resample to subject <qtarget>, e.g. fsaverage, and smooth at differnt levels)
To manually specify the subject specific output in <template>/surf/<hemi>.<name> :
--name-avg <name> In/out filename (without hemi ?h) for temporal average (default: long.<meas>-avg.mgh)
--name-rate <name> In/out filename for the rate (diff per time) (default: long.<meas>-rate.mgh)
--name-pc1fit <name> In/out filename for percent change (w.r.t. time 1 fit) (default: long.<meas>-pc1fit.mgh)
--name-pc1 <name> In/out filename for percent change (w.r.t. time 1) (default: long.<meas>-pc1.mgh)
--name-spc <name> In/out filename for sym. pct. change (w.r.t. average) (default: long.<meas>-spc.mgh)
--out-stack <name> Store stacked measure file <template>/surf/<hemi>.<name> (default: long.<meas>-stack.mgh)
--out-label <name> Store intersected cortex in <template>/label/<hemi>.<name>.label (default: long.cortex)
To ouput stacked/intersected results on the qtarget (all require --qcache) :
--isec-labels <name> Intersect labels on <qtarget> (usually cortex labels)
--stack-avg <name> Output stacked avg maps on <qtarget> for all fwhm levels
--stack-rate <name> Output stacked rate maps on <qtarget> for all fwhm levels
--stack-pc1fit <name> Output stacked pc1fit maps on <qtarget> for all fwhm levels
--stack-pc1 <name> Output stacked pc1 maps on <qtarget> for all fwhm levels
--stack-spc <name> Output stacked spc maps on <qtarget> for all fwhm levels
"""
HELPTEXT = """
SUMMARY
=======
Computes slope maps (e.g. of thickness) in a longitudinal study.
The slope is computed within subject from the longitudinally processed
results (taken from the <tpNid>.long.<template> directories) and the
output is written into the subjects <template>/surf directory for further
processing (e.g. group analysis). Will also output the intersections of
all within-subject cortex labels in <template>/label/?h.long.cortex.label
if an input label is not specified with --inlabel <name> or if no different
output name is specified with --outlabel <name> .
DETAILS
=======
QDEC.TABLE
==========
Pass a qdec table file, where the first 2 columns need to be 'fsid fsid-base'.
fsid is the id of the individual time points an 'fsid-base' the template/base
id (grouping the timepoints that belong to the same subject). By default the
third column is taken as the time variable, but this can be overwritten with
--time <name>.
QDEC.TABLE-EXAMPLE:
fsid fsid-base age weight IQ
Elmo_1 Elmo 3 10 1000
#Elmo_2 Elmo 3.5 15 1100
Elmo_3 Elmo 4 20 1300
Snuffy_1 Snuffy 20 40 1100
Snuffy_2 Snuffy 21 45 1200
Bert_1 Bert 8 25 2000
Bert_2 Bert 9 30 2500
Bert_3 Bert 9.9 34 2400
If you like, you can create an independent qdec table file for each subject,
containing only the time points of that subject (e.g. to process stuff
simultaneously). You can also run left and right hemisphere in parallel.
ANALYSIS
========
Choose one or more of the following options:
* 'do-rate' for rate, this will compute the slope of a linear fit across
time within each subject. Depending on the time variable this will
yield the thinning in mm/time (if the time variable is measured in
years, such as age: mm/year).
* 'do-pc1' is the percent change, this is the rate normalized by the
measure at the first time point times 100, e.g. percent thinning per
year.
* 'do-pc1fit' is the percent change with respect to the value at
tp1, evalutated from the linear fit. This is the rate normalized by the
estimated measure at the first time point times 100, e.g. percent
thinning per year. Using the linear fit instead of the measure at time
point 1 directly reduces noise as tp1 data is more noisy than the fit,
therefore we recommend pc1fit over pc1.
* 'do-spc' for symmetrized percent change. Here we normalize by the temporal
average instead of taking it from the first time point. The average
is computed from the linear fit at the middle of the time interval.
This is a symmetric 'percent thinning per year' and more robust, as
TP1 can be an outlier.
* 'do-avg' for output of the temporal average (linear fit evaluated at mid
time).
NOTE that for percent change (pc1fit, pc1 and spc) prior smoothing (--fwhm)
is highly recommended as the division can produce very large values!
WITHIN SUBJECT OUTPUT
=====================
The output of the within subject surface maps will be written into the
template's surf directory:
<template>/surf/<hemi>.<NAME>(.fwhm<fwhm>).mgh
For each of the above analysis, the <NAME> can be manually specified
using --name-rate, --name-pc1fit, --name-pc1, --name-spc, --name-avg.
Defaults are: long.<meas>-rate, long.<meas>-pc1 etc.
It is possible to also store the stack of the measure across time in each
subject using -do-stack. The default is destination is
<template>/surf/<hemi>.<NAME>(.fwhm<fwhm>).mgh
with <NAME> = 'long.<meas>-stack' which can be overwritten with
--out-stack <NAME>
CORTEX LABELS
=============
Create within subject intersections:
It is recommended to specify --do-label which will intersect the cortex label
of all time points within each subject. The intersected label is saved to
<template>/label/<HEMI>.<OUT_LABEL> (default: <OUT_LABEL>=long.cortex).
The <OUT_LABEL> part can be manually specified by --out-label <OUT_LABEL>.
Map to common target subject (qcache, see also below):
If additionally --qcache and --isec-labels <name> is specified, these labels
will be mapped and intersected on the <qtarget> and output will be written
to <name> (where <name> is the full path and filename of the output label).
Use existing labels:
An existing label can be used and is specified with --in-label <name> where
<name> is the label name without the hemisphere as found in each subjects
template label dir: template/label/<hemi>.<name>
for example 'long.cortex.label' to use previously created intersected labels
or just 'cortex.label' to use the label from each template.
Only intersect existing qtarget labels:
If labels were mapped to qtarget (with --qcache) before and need to be
intersected again (e.g. if subjects were removed) use
--in-label <name>.<qtarget>.label --qcache <qtarget> --iscec-labels <outlabel>
but don't specify any do-??? flags as within subject analysis will fail
using the qtarget labels.
QCACHE
======
In order to run a group analysis with QDEC all surface measures need to be
resampled to some target. If you plant to use QDEC, specify --qcache <qtarget>.
You can select the qtarget (needs to be in the subjectsdir), usually fsaverage
is fine, which is automatically created in the subjectsdir during recon-all.
If --qcache <qtarget> is specified, it will resample all the analyses
specified with --do-??? (and the labels if --do-label) to the qtarget.
It will also work on different smoothing levels if --fwhm is NOT specified:
0,5,10,15,20,25
If --fwhm is specified it will only work on the specified level.
Smoothing will be done in the subject space before the division (spc or pc1).
Output can be found in <template>/surf folder under each subject template.
<template>/surf/<HEMI>.<NAME>(.fwhm<fwhm>).mgh (for subject space)
<template>/surf/<HEMI>.<NAME>(.fwhm<fwhm>).<QTARGET>.mgh
If cortex labels are intersected (--do-label), it will also be mapped to
the qtarget and output as:
<template>/label/<HEMI>.<NAME>(.fwhm<fwhm>).<QTARGET>.mgh
Stacking all maps on QTARGET:
You can specify --stack-??? to stack all maps on qtarget, this is needed
if you plan to run your own GLM's (mri_glmfit). It is possible to stack
previously created qtarget maps, simply do not run any subject analysis
with do-???.
EXAMPLES
========
Simply create avg, rate, pc1, spc for each subject. Also will output
the stacked thickness files and will intersect the cortex label. Here we
work on the left hemisphere (lh) and use smoothing of 15 fwhm. Note that
'age' is the column secifying the age in the long.qdec.table.dat file:
long_mris_slopes --qdec ./long.qdec.table.dat --meas thickness --hemi lh \
--do-avg --do-rate --do-pc1fit --do-spc --do-stack --do-label \
--fwhm 15 --time age --sd ./
To create Qdec cache on fsaverage (several fwhm levels):
long_mris_slopes --qdec ./long.qdec.table.dat --meas thickness --hemi lh \
--do-avg --do-rate --do-pc1fit --do-spc --do-stack --do-label \
--time age --qcache fsaverage --sd ./
To stack results for spc and fwhm 15:
long_mris_slopes --qdec ./long.qdec.table.dat --meas thickness --hemi lh \
--fwhm 15 --qcache fsaverage --stack-spc ./lh.thickness-spc.stack.mgh \
--sd ./
To intersect previously mapped cortex labels:
long_mris_slopes --qdec ./long.qdec.table.dat --hemi lh \
--qcache fsaverage --in-label long.cortex.fsaverage.label \
--isec-labels ./lh.long.isec.fsaverage.cotex.label --sd ./
REFERENCES
==========
Highly Accurate Inverse Consistent Registration: A Robust Approach,
M. Reuter, H.D. Rosas, B. Fischl. NeuroImage 53(4), 1181-1196, 2010.
http://dx.doi.org/10.1016/j.neuroimage.2010.07.020
http://reuter.mit.edu/papers/reuter-robreg10.pdf
Avoiding Asymmetry-Induced Bias in Longitudinal Image Processing,
M. Reuter, B. Fischl. NeuroImage 57(1), 19-21, 2011.
http://dx.doi.org/10.1016/j.neuroimage.2011.02.076
http://reuter.mit.edu/papers/reuter-bias11.pdf
Within-Subject Template Estimation for Unbiased Longitudinal Image Analysis.
M. Reuter, N.J. Schmansky, H.D. Rosas, B. Fischl.
NeuroImage 61(4), 1402-1418, 2012.
http://dx.doi.org/10.1016/j.neuroimage.2012.02.084
http://reuter.mit.edu/papers/reuter-long12.pdf
"""
def options_parse():
"""
Command Line Options Parser for long_mris_slopes
initiate the option parser and return the parsed object
"""
parser = optparse.OptionParser(version='$Id: long_mris_slopes,v 1.38 2012/05/30 22:57:25 mreuter Exp $', usage=HELPTEXT)
# help text
# h_fsgd = '(REQUIRED) fsdg file specifying the subjects and time points'
h_qdec = '(REQUIRED) qdec table file specifying the subjects and time points'
h_meas = '(REQUIRED) the surface input measure (e.g. thickness)'
h_hemi = '(REQUIRED) run one hemisphere: lh or rh or both'
h_sd = '(REQUIRED) full path to FreeSurfer subjects dir'
h_doavg = 'compute and output the temporal average (recommended)'
h_dorate = 'compute and output the rate (recommended)'
h_dopc1fit = 'compute and output the pct. change (w.r.t. tp1 from linear fit)'
h_dopc1 = 'compute and output the pct. change (w.r.t. tp1)'
h_dospc = 'compute and output the sym. pct. change (w.r.t. temp. average) (recommended)'
h_dostack = 'save the stacked within subject file (time series)'
h_dolabel = 'compute and output intersected cortex label (recommended)'
h_resid = 'residual tp (pass 1 for tp1) to export (default no export)'
h_name_avg = 'filename (without hemi ?h) to store temporal average in <template>/surf/<hemi>.<NAME_AVG>.mgh (default: long.<meas>-avg)'
h_name_rate = 'filename (without hemi ?h) to store rate in <template>/surf/<hemi>.<NAME_RATE>.mgh (default: long.<meas>-rate)'
h_name_pc1fit = 'filename (without hemi ?h) to store pct. change (to TP1fit) in <template>/surf/<hemi>.<NAME_PC1FIT>.mgh (default: long.<meas>-pc1fit)'
h_name_pc1 = 'filename (without hemi ?h) to store pct. change (to TP1) in <template>/surf/<hemi>.<NAME_PC1>.mgh (default: long.<meas>-pc1)'
h_name_spc = 'filename (without hemi ?h) to store sym. pct. change in <template>/surf/<hemi>.<NAME_SPC>.mgh (default: long.<meas>-spc)'
h_name_resid = 'filename (without hemi ?h) to store residual in <template>/surf/<hemi>.<NAME_RESID>.mgh (default: long.<meas>-resid<resid>), requires --resid <int> '
h_outstack = 'filename to store stacked measure file <template>/surf/<hemi>.<OUT_STACK>.mgh (default: long.<meas>-stack)'
h_outlabel = 'filename to store within-subject intersected cortex labels in <template>/label/<hemi>.<OUT_LABEL>.label (default: long.cortex)'
h_qcache = 'create cache for qdec (resample to subject <QCACHE>, e.g. fsaverage, and smooth at differnt levels. Also specify qcache output names (see below))'
h_fwhm = 'smoothing at specific fwhm (required for pct. change)'
h_nosmooth = 'do not smooth the data (overwrite requirement to not smooth pct change computations)'
h_time = 'variable name for time column variable (e.g. age) in qdec table'
h_generictime = 'time points are ordered in qdec file, assume time=1,2,3...'
h_inlabel = 'use pre-existing label for smoothing and to mask the output'
h_jac = 'use this flag when mapping area or volume maps to correct Jacobian (passed to mri_surf2surf)'
h_iseclabels= 'intersect labels on <qtarget> (usually cortex labels), requires --qcache'
h_stackavg = 'output stacked avg maps on <qtarget> for all fwhm levels, requires --qcache'
h_stackrate = 'output stacked rate maps on <qtarget> for all fwhm levels, requires --qcache'
h_stackpc1fit = 'output stacked pc1fit maps on <qtarget> for all fwhm levels, requires --qcache'
h_stackpc1 = 'output stacked pc1 maps on <qtarget> for all fwhm levels, requires --qcache'
h_stackspc = 'output stacked spc maps on <qtarget> for all fwhm levels, requires --qcache'
h_stack_resid = 'output stacked residual maps on <qtarget> for all fwhm levels, requires --qcache'
# h_outformat = 'Output format can be either mgh (default) or curv'
# Add options
# parser.add_option('--fsgd', dest='fsgd', help=h_fsgd)
# Sepcify inputs
group = optparse.OptionGroup(parser, "Required Arguments")
group.add_option('--qdec', dest='qdec', help=h_qdec)
group.add_option('--meas', dest='meas', help=h_meas)
group.add_option('--hemi', dest='hemi', help=h_hemi, choices=('lh','rh','both'))
group.add_option('--sd', dest='sd' , help=h_sd)
parser.add_option_group(group)
# do computations:
group = optparse.OptionGroup(parser, "Computations", "Select one or more of the following options:")
group.add_option('--do-avg' , action='store_true', dest='do_avg' , help=h_doavg , default=False)
group.add_option('--do-rate' , action='store_true', dest='do_rate' , help=h_dorate , default=False)
group.add_option('--do-pc1fit' , action='store_true', dest='do_pc1fit' , help=h_dopc1fit , default=False)
group.add_option('--do-pc1' , action='store_true', dest='do_pc1' , help=h_dopc1 , default=False)
group.add_option('--do-spc' , action='store_true', dest='do_spc' , help=h_dospc , default=False)
group.add_option('--do-stack', action='store_true', dest='do_stack', help=h_dostack, default=False)
group.add_option('--resid', dest='resid', type="int", help=h_resid, default=-1)
group.add_option('--do-label', action='store_true', dest='do_label', help=h_dolabel, default=False)
group.add_option('--qcache' , dest='qcache' , help=h_qcache)
parser.add_option_group(group)
# parameters:
group = optparse.OptionGroup(parser, "Parameters")
group.add_option('--fwhm', dest='fwhm', help=h_fwhm)
group.add_option('--nosmooth', action='store_true', dest='nosmooth', help=h_nosmooth, default=False)
group.add_option('--time', dest='time', help=h_time)
group.add_option('--generic-time', action='store_true', dest='generic_time', help=h_generictime, default=False)
group.add_option('--in-label' , dest='in_label' , help=h_inlabel)
group.add_option('--jac', action='store_true', dest='jac', help=h_jac, default=False)
parser.add_option_group(group)
# overwrite default names:
group = optparse.OptionGroup(parser, "Within-Subject Output","Select these to overwrite the default names:")
group.add_option('--name-avg' , dest='name_avg' , help=h_name_avg)
group.add_option('--name-rate' , dest='name_rate' , help=h_name_rate)
group.add_option('--name-pc1fit' , dest='name_pc1fit' , help=h_name_pc1fit)
group.add_option('--name-pc1' , dest='name_pc1' , help=h_name_pc1)
group.add_option('--name-spc' , dest='name_spc' , help=h_name_spc)
group.add_option('--name-resid' , dest='name_resid' , help=h_name_resid)
group.add_option('--out-stack' , dest='out_stack', help=h_outstack)
group.add_option('--out-label' , dest='out_label', help=h_outlabel, default="long.cortex")
parser.add_option_group(group)
# stuff on common qcache target:
group = optparse.OptionGroup(parser, "Outputs on Common <qcache> Target","To specify outputs when running --qcache <id>")
group.add_option('--isec-labels', dest='isec_labels', help=h_iseclabels)
group.add_option('--stack-avg' , dest='stack_avg' , help=h_stackavg)
group.add_option('--stack-rate' , dest='stack_rate' , help=h_stackrate)
group.add_option('--stack-pc1fit' , dest='stack_pc1fit' , help=h_stackpc1fit)
group.add_option('--stack-pc1' , dest='stack_pc1' , help=h_stackpc1)
group.add_option('--stack-spc' , dest='stack_spc' , help=h_stackspc)
group.add_option('--stack-resid' , dest='stack_resid' , help=h_stack_resid)
parser.add_option_group(group)
(options, args) = parser.parse_args()
# extensive error checks
# if options.fsgd is None and
if options.qdec is None:
# print 'ERROR: Specify --fsgd or --qedc'
parser.print_help()
print '\nERROR: Specify --qedc (and other required arguments)\n'
sys.exit(1)
# if options.fsgd is not None and options.qdec is not None:
# print 'ERROR: Specify either --fsgd or --qedc, not both'
# sys.exit(1)
if options.hemi is None:
print 'ERROR: Specify --hemi (e.g. \'lh\')\n'
sys.exit(1)
if options.do_avg and options.name_avg is None:
options.name_avg = 'long.'+options.meas+'-avg'
# elif options.name_avg is not None:
# options.do_avg = True
if options.do_rate and options.name_rate is None:
options.name_rate = 'long.'+options.meas+'-rate'
# elif options.name_rate is not None:
# options.do_rate = True
if options.do_pc1fit and options.name_pc1fit is None:
options.name_pc1fit = 'long.'+options.meas+'-pc1fit'
if options.do_pc1 and options.name_pc1 is None:
options.name_pc1 = 'long.'+options.meas+'-pc1'
# elif options.name_pc1 is not None:
# options.do_pc1 = True
if options.do_spc and options.name_spc is None:
options.name_spc = 'long.'+options.meas+'-spc'
# elif options.name_spc is not None:
# options.do_spc = True
if options.do_stack and options.out_stack is None:
options.out_stack = 'long.'+options.meas+'-stack'
# elif options.out_stack is not None:
# options.do_stack = True
if options.resid > 0 and options.name_resid is None:
options.name_resid = 'long.'+options.meas+'-resid'+str(options.resid)
if options.resid < 1 and options.name_resid is not None:
print 'ERROR: Please also specify the tp number to export residuals (starting at 1) with --resid !\n'
sys.exit(1)
do_analysis = options.do_avg or options.do_pc1fit or options.do_pc1 or options.do_spc or options.do_rate or (options.resid > 0)
if options.meas is None and do_analysis:
print 'ERROR: Specify --meas (e.g. \'thickness\')\n'
sys.exit(1)
if options.sd is None:
print 'ERROR: Specify the subject dir with --sd <fullpath>\n'
sys.exit(1)
if options.meas is None: options.meas = "label" #needed later to create tmp dir
do_stacking = options.stack_spc is not None or options.stack_pc1fit is not None or options.stack_pc1 is not None or options.stack_rate is not None or options.stack_avg is not None
if do_stacking and options.qcache is None:
print 'ERROR: --stack-??? requires target subject (e.g. fsaverage) as specified with --qcache\n'
sys.exit(1)
# if we run analysis, but no label is created or passed
if do_analysis and not options.do_label and options.in_label is None:
print 'ERROR: either create intersected labels with --do-label or specify input label names with --in-label\n'
sys.exit(1)
# if labels should be created and are also passed:
if options.do_label and options.in_label is not None:
print 'ERROR: either create intersected labels with --do-label or specify input label names with --in-label\n'
# switch on do_label if output is specified but no input
if options.isec_labels is not None and options.in_label is None:
options.do_label = True
if not do_analysis and not options.do_label and not do_stacking and options.isec_labels is None:
print 'ERROR: Analysis type should be specified, use one or more of --do-avg, --do-rate, --do-pc1, --do-spc, --do-label, --stack-??? or --isec-labels\n'
sys.exit(1)
if options.fwhm is not None and options.nosmooth:
print 'ERROR: either specify --fwhm <int> or --nosmooth\n'
sys.exit(1)
if (options.name_pc1fit is not None or options.name_pc1 is not None or options.name_spc is not None) and options.fwhm is None and not options.nosmooth and options.qcache is None:
print 'ERROR: when computing SPC or PC1 (pct. change) you should specify smoothing --fwhm <int>. If you really do not want to smooth pass --nosmooth\n'
sys.exit(1)
if options.isec_labels is not None and options.qcache is None:
print 'ERROR: --isec-labels requires target subject (e.g. fsaverage) as specified with --qcache\n'
sys.exit(1)
return options
def write_fsgd(filename,subjects_tp_map,subject,timeidx):
"""
Write the fsgd to disk. Initialize the writer class.
"""
fp = open(filename, 'w')
fp.write('GroupDescriptorFile 1\n')
fp.write('Title MyLongStudy\n')
fp.write('Class '+subject)
fp.write('Variables Time')
for tp in subjects_tp_map[subject]:
fp.write('Input '+tp[0]+' '+subject+' '+tp[timeidx]+'\n')
fp.write('\n')
fp.close()
def run_cmd(cmd,err_msg):
"""
execute the comand
"""
print cmd+'\n'
args = shlex.split(cmd)
retcode = subprocess.call(args)
if retcode != 0 :
print 'ERROR '+str(retcode)+' : '+err_msg
sys.exit(1)
print '\n'
def intersect_labels(labellist,outlabel):
# will intersect the labels in labelist (full paths)
# and write the output to outlabel (full path)
if os.path.exists(outlabel):
os.remove(outlabel)
first = True
for label in labellist:
if not os.path.exists(label):
print 'ERROR: '+str(label)+' does not exist!'
sys.exit(1)
if first:
cmd = 'cp '+label+' '+outlabel
first = False
else:
cmd = 'mris_label_calc intersect '+label+' '+outlabel+' '+outlabel
run_cmd(cmd,'mris_label_calc intersect did not work?')
if not os.path.exists(outlabel):
print 'ERROR: Cannot write '+str(outlabel)+' (write permissions?)!'
sys.exit(1)
def stack_map(srcid,trgid,src_type,trg_type,inlist,outlist,tmpdir,hemi,jac):
if len(inlist) != len(outlist):
print 'ERROR: stack_map: in and outlist have different lengths!'
sys.exit(1)
if len(inlist) ==0:
print 'ERROR: stack_map: inlist should have at least one element!'
sys.exit(1)
# jac is for mapping area and volume maps
jacflag = ''
if jac:
jacflag = ' --jac'
if len(inlist) == 1:
# map stack to target
trg1=os.path.join(tmpdir,hemi+'.tmp0000.mgh')
cmd = 'mri_surf2surf'+jacflag+' --hemi '+hemi+' --srcsubject '+srcid+' --sval '+inlist[0]+' --src_type mgh --trgsubject '+trgid+' --tval '+trg1+' --trg_type mgh'
run_cmd(cmd,'mri_surf2surf mapping single file to qcache target did not work?')
else:
allin = " ".join(inlist)
# stack all measure files
allstack=os.path.join(tmpdir,'allstack.mgh')
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking all outputs')
# map stack to target
trgstack=os.path.join(tmpdir,'trgstack.mgh')
cmd = 'mri_surf2surf'+jacflag+' --hemi '+hemi+' --srcsubject '+srcid+' --sval '+allstack+' --src_type mgh --trgsubject '+trgid+' --tval '+trgstack+' --trg_type mgh'
run_cmd(cmd,'mri_surf2surf mapping stack to qcache target did not work?')
# split the stack
outframe=os.path.join(tmpdir,hemi+'.tmp.mgh')
cmd = 'mri_convert '+trgstack+' '+outframe+' --split'
run_cmd(cmd,'mri_convert splitting qcache stack did not work?')
# mv elements to final destination
i=0
for fn in outlist:
istr= '%04d' % i
tmpfn=os.path.join(tmpdir,hemi+'.tmp'+istr+'.mgh')
if not os.path.exists(tmpfn):
print 'ERROR: '+str(tmpfn)+' does not exist!'
sys.exit(1)
cmd = 'mv '+tmpfn+' '+fn
run_cmd(cmd,'mv '+tmpfn+' '+fn+' did not work?')
i=i+1
if __name__=="__main__":
# Command Line options and error checking done here
options = options_parse()
slopelogger.debug('-- The options you entered --')
slopelogger.debug(options)
subjectsdir = ''
# Parse the stats files
# if options.fsgd is not None:
# print 'Parsing the fsgd file: '+options.fsgd
# try:
# slopelogger.debug('Processing file ' + options.fsgd)
# fsgdparse = FsgdParser(options.fsgd)
# subjects_tp_map, variables, defaultvar = fsgdparse.parse()
# except BadFileError, e:
# print 'ERROR: fsgd file '+str(e)+' not found!'
# sys.exit(1)
# else:
print 'Parsing the qdec table: '+options.qdec
try:
slopelogger.debug('Processing file ' + options.qdec)
qdectable = LongQdecTable(options.qdec)
#subjects_tp_map, variables, subjectdir = qdecparse.parse()
except BadFileError, e:
print 'ERROR: qdec table '+options.qdec+' not found or empty?'
sys.exit(1)
# make sure we have a long table containing the bases
if qdectable.cross:
print '\nERROR: qdec table '+options.qdec+' is cross sectional\n (2nd column not \'fsid-base\')!\n'
sys.exit(1)
# use the first column by default for time variable
varidx = 1
# get other variables:
variables=qdectable.variables
if len(variables) < 1:
varidx = -1
if not options.generic_time:
print '\nERROR: qdec table '+options.qdec+' needs 3rd column with time value,'
print ' e.g. age or time since baseline scan ...!\n'
print ' Or pass --generic-time assuming tps are ordered 1,2,3...\n'
sys.exit(1)
# if time variable is passed on command line, make sure it is part of variables
if not options.time is None:
defaultvar = options.time
# compute correct index (starting with 1, 0 is the tpID)
foundidx = False
for index in (i for i in xrange(len(variables)) if variables[i].upper()==defaultvar.upper()):
varidx = index
foundidx = True
#print 'found: '+str(varidx)+' '+variables[varidx]
break
if not foundidx:
print '\nERROR: DefaultVariable \''+str(defaultvar)+'\' not found in variables: '+str(variables)+'\n'
sys.exit(1)
varidx = varidx +1;
# maybe later check if time column is really a float?
# if env is set, overwrite info from file (if it was passed in qdec)
subjectsdir = options.sd
if subjectsdir is None:
subjectsdir = os.getenv('SUBJECTS_DIR')
if subjectsdir is None:
subjectsdir = qdectable.subjectsdir
if subjectsdir is None:
print'\nERROR: no subjects dir specified, use --sd <fullpath>'
sys.exit(1)
print '\nWorking in SUBJECTS_DIR: '+subjectsdir+'\n'
os.environ['SUBJECTS_DIR'] = subjectsdir
if options.qcache is not None:
# check if qcache target exists
qsdir = os.path.join(subjectsdir,options.qcache)
if not os.path.exists(qsdir):
print 'ERROR: Qcache Target dir '+str(qsdir)+' not found!'
sys.exit(1)
hemis = [ options.hemi ]
if hemis[0] == "both": hemis = [ "lh" , "rh" ]
for hemi in hemis:
print
print 'Working on hemi: '+hemi
print
# process
retcode = 0
tlabels = []
tavg = {}
trate = {}
tpc1 = {}
tpc1fit = {}
tspc = {}
tresid = {}
for subjectid, tplist in qdectable.subjects_tp_map.items():
print '\nSubject-Template: '+subjectid
# check if basedir exists
basedir = os.path.join(subjectsdir,subjectid)
if not os.path.exists(basedir):
print '\nERROR: Template dir '+str(basedir)+' does not exist!\n '
sys.exit(1)
basesurfdir=os.path.join(basedir,'surf')
#if options.qcache is not None:
# basesurfdir=os.path.join(basedir,'surf','qcache')
# if not os.path.exists(basesurfdir):
# os.mkdir(basesurfdir)
# extract ids and age data:
tpids = [entry[0]+'.long.'+subjectid for entry in tplist]
if varidx > 0:
times = [float(entry[varidx]) for entry in tplist]
else:
times = [float(x) for x in range(1,len(tplist)+1)]
num = len(times)
meant = sum(times) / num
print '\n\nINFO: '+str(num)+' TPs in '+subjectid+' , mean age: '+str(meant)+'\n'
#print str(times)
#sys.exit(0)
# intersect cortex ?
label=""
if options.do_label:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" Intersecting Within-Subject Cortex Label\n"
label = hemi+'.'+options.out_label+'.label'
label = os.path.join(basedir,'label',label)
llist = [ os.path.join(subjectsdir,entry,'label',hemi+'.cortex.label') for entry in tpids]
# ouput written to base: /label/<hemi>.<outlabel>.label , default outlabel=long.cortex
intersect_labels(llist,label)
# map label of this subject to common qcache target:
if options.qcache is not None:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" mapping label to QCache\n"
label = hemi+'.'+options.out_label+'.label'
label = os.path.join(basedir,'label',label)
outlabel = hemi+'.'+options.out_label+'.'+options.qcache+'.label'
outlabel = os.path.join(basedir,'label',outlabel)
cmd = 'mri_label2label --srclabel '+label+' --srcsubject '+subjectid+' --trglabel '+outlabel+' --trgsubject '+options.qcache
cmd = cmd+' --regmethod surface --hemi '+hemi
run_cmd(cmd,'mri_label2label mapping the label failed?')
tlabels.append(outlabel)
if options.in_label is not None:
label = os.path.join(basedir,'label',options.in_label)
if not os.path.exists(label):
label = hemi+'.'+options.in_label
label = os.path.join(basedir,'label',label)
if not os.path.exists(label):
label = hemi+'.'+options.in_label+'.label'
label = os.path.join(basedir,'label',label)
if not os.path.exists(label):
print '\nERROR: cannot find '+options.in_label+' or '+label+'\n'
sys.exit(1)
tlabels.append(label)
# setup fwhms (either specified or full list):
fwhms = [ 0 ]
if options.qcache is not None and options.fwhm is None:
fwhms = [ 0 , 5, 10, 15, 20, 25 ]
elif options.fwhm is not None:
fwhms = [ options.fwhm ]
# create lists of ouput filenames for stacking on qcache (depending on type)
if options.qcache is not None:
for fwhm in fwhms:
fwhmstr = '.fwhm'+str(fwhm)
if options.stack_resid is not None:
outtname= os.path.join(basesurfdir,hemi+'.'+options.name_resid+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in tresid:
tresid[fwhm] = []
tresid[fwhm].append(outtname)
if options.stack_spc is not None:
outtname= os.path.join(basesurfdir,hemi+'.'+options.name_spc+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in tspc:
tspc[fwhm] = []
tspc[fwhm].append(outtname)
if options.stack_pc1fit is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_pc1fit+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in tpc1fit:
tpc1fit[fwhm] = []
tpc1fit[fwhm].append(outtname)
if options.stack_pc1 is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_pc1+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in tpc1:
tpc1[fwhm] = []
tpc1[fwhm].append(outtname)
if options.stack_rate is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_rate+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in trate:
trate[fwhm] = []
trate[fwhm].append(outtname)
if options.stack_avg is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_avg+fwhmstr+'.'+options.qcache+'.mgh')
if fwhm not in tavg:
tavg[fwhm] = []
tavg[fwhm].append(outtname)
# if no subject base processing, continue:
if not (options.do_pc1fit or options.do_pc1 or options.do_avg or options.do_spc or options.do_rate or options.resid > 0):
continue #skip within subject processing and mapping to target
# check if 2 or more time points
if len(tplist) < 2 :
print 'ERROR: '+str(basedir)+' must have at least 2 time points!'
sys.exit(1)
# list of thickness files:
measures=[os.path.join(subjectsdir,entry,'surf',hemi+'.'+options.meas) for entry in tpids]
all = " ".join(measures)
# create tmpdir:
prefix = './tmp-'+subjectid+'_'+hemi+'_'+options.meas+'_'
dirname = tempfile.mkdtemp('',prefix,'')
if not os.path.exists(dirname):
print 'ERROR: tmp dir '+str(dirname)+' cannot be created (write permissions?)!'
sys.exit(1)
# stack measure maps:
meas_target = os.path.join(dirname,hemi+'.long.'+options.meas+'-stack.mgh')
if options.do_stack:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" Stacking Within-Subject Maps\n"
meas_target = os.path.join(basesurfdir,hemi+'.'+options.out_stack)
if meas_target[-4:] != '.mgh':
meas_target=meas_target+'.mgh'
cmd = 'mri_concat '+all+' --o '+meas_target
run_cmd(cmd,'mri_concat stacking did not work?')
# # 2 tps, direct computation
# if num == 2 :
#
# # sort
# if times[0] > times[1] :
# times[0] , times[1] = times[1] , times[0]
# measures[0] , measures[1] = measures[1], measures[0]
# tpids[0], tpids[1] = tpids[1], tpids[0]
#
# labelpara=""
# if label != "" :
# labelpara = ' --label '+label
#
# # create difference:
# beta1 = os.path.join(dirname,hemi+'.rate.'+options.meas)
# cmd = 'mris_calc -o '+beta1+labelpara+' '+measures[1]+' sub '+measures[0]
# run_cmd(cmd,'mri_calc difference did not work?')
# cmd = 'mris_calc -o '+beta1+labelpara+' '+beta1+' div '+str(times[1]-times[0])
# run_cmd(cmd,'mri_calc div did not work?')
#
# # create average for SPC or avg
# beta0 = os.path.join(dirname,hemi+'.avg.'+options.meas)
# cmd = 'mris_calc -o '+beta0+labelpara+' '+measures[0]+' add '+measures[1]
# run_cmd(cmd,'mri_calc sum did not work?')
# cmd = 'mris_calc -o '+beta0+' '+beta0+' mul 0.5'
# run_cmd(cmd,'mri_calc mul 0.5 did not work?')
#
# # convert to mgh format (to imitate glm)
# beta0old = beta0
# beta0 = os.path.join(dirname,'beta0.mgh')
# beta1old = beta1
# beta1 = os.path.join(dirname,'beta1.mgh')
# cmd = 'mri_convert '+beta0old+' '+beta0
# run_cmd(cmd,'mri_convert did not work?')
# cmd = 'mri_convert '+beta1old+' '+beta1
# run_cmd(cmd,'mri_convert did not work?')
#
#
# else: # more than 2 time points, run GLM
if 1==1:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" Running Within-Subject GLM\n"
#write X-matrix (times):
x_target = os.path.join(dirname,'X-long.mat')
print 'Writing '+x_target+' ...\n'
if os.path.exists(x_target):
os.remove(x_target)
xfp = open(x_target, 'w')
for time in times:
xfp.write('1 '+str(time-times[0])+'\n')
#xfp.write('1 '+str(time-meant)+'\n')
xfp.close()
# run glm in tmp dir:
zerodof=""
if num==2:
zerodof=" --allow-zero-dof "
saveeres=""
if options.resid > 0:
saveeres=" --eres-save "
labelpara=""
if label != "" :
labelpara = ' --label '+label
glmdir=os.path.join(dirname,'glm')
cmd = 'mri_glmfit --y '+meas_target+' --X '+x_target+zerodof+saveeres+' --no-contrasts-ok --surf '+subjectid+' '+hemi+labelpara+' --glmdir '+glmdir
run_cmd(cmd,'mri_glmfit did not work?')
# harvest results (in beta.mgh)
betafn = os.path.join(glmdir,'beta.mgh' )
beta0 = os.path.join(dirname,'beta0.mgh')
beta1 = os.path.join(dirname,'beta1.mgh')
betam = os.path.join(dirname,'betam.mgh')
outresid = ""
if not os.path.exists(betafn):
print 'ERROR: GLM results '+str(betafn)+' does not exist!'
sys.exit(1)
# split beta
cmd = 'mri_convert --frame 0 '+betafn+' '+beta0
run_cmd(cmd,'mri_convert split frames 0 did not work?')
cmd = 'mri_convert --frame 1 '+betafn+' '+beta1
run_cmd(cmd,'mri_convert split frames 1 did not work?')
# evaluate fit at mean age: betam = beta0 + (meanage-baseage) beta1
cmd = 'mris_calc -o '+betam+' '+beta1+' mul '+str(meant-times[0])
run_cmd(cmd,'mris_calc compute betam (mul) problem?')
cmd = 'mris_calc -o '+betam+' '+betam+' add '+beta0
run_cmd(cmd,'mris_calc compute betam (add) problem?')
# get residuals if needed
if options.resid > 0:
eres = os.path.join(glmdir,'eres.mgh' )
outresid = os.path.join(dirname,options.name_resid+'.mgh')
cmd = 'mri_convert --frame '+str(options.resid - 1)+' '+eres+' '+outresid
run_cmd(cmd,'mri_convert split residual frames '+str(options.resid)+' did not work?')
# either of the cases:
allout=[]
alltout=[]
for fwhm in fwhms:
fwhmstr = ''
if options.qcache is not None or fwhm > 0 :
fwhmstr = '.fwhm'+str(fwhm)
# do smoothing
beta0fwhm = ''
beta1fwhm = ''
betamfwhm = ''
tp1fwhm = ''
residfwhm = ''
# Smooth data on subject level (before computations)
if fwhm > 0:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" Smoothing at "+str(fwhm)+"\n"
smooth = ' --fwhm '+str(fwhm)+' --smooth-only'
labelpara=""
if label != "" :
labelpara = ' --label '+label
beta0fwhm = os.path.join(dirname,'beta0'+fwhmstr+'.mgh')
cmd = 'mris_fwhm --hemi '+hemi+' --s '+subjectid+' --i '+beta0+labelpara+smooth+' --o '+beta0fwhm
run_cmd(cmd,'mris_fwhm smoothing beta0 did not work?')
beta1fwhm = os.path.join(dirname,'beta1'+fwhmstr+'.mgh')
cmd = 'mris_fwhm --hemi '+hemi+' --s '+subjectid+' --i '+beta1+labelpara+smooth+' --o '+beta1fwhm
run_cmd(cmd,'mris_fwhm smoothing beta1 did not work?')
betamfwhm = os.path.join(dirname,'betam'+fwhmstr+'.mgh')
cmd = 'mris_fwhm --hemi '+hemi+' --s '+subjectid+' --i '+betam+labelpara+smooth+' --o '+betamfwhm
run_cmd(cmd,'mris_fwhm smoothing betam did not work?')
if options.do_pc1:
tp1fwhm = os.path.join(dirname,'tp1'+fwhmstr+'.mgh')
cmd = 'mris_fwhm --hemi '+hemi+' --s '+tpids[0]+' --i '+measures[0]+labelpara+smooth+' --o '+tp1fwhm
run_cmd(cmd,'mris_fwhm smoothing tp1 did not work?')
if options.resid > 0:
residfwhm = os.path.join(dirname,'resid'+fwhmstr+'.mgh')
cmd = 'mris_fwhm --hemi '+hemi+' --s '+subjectid+' --i '+outresid+labelpara+smooth+' --o '+residfwhm
run_cmd(cmd,'mris_fwhm smoothing resid did not work?')
else:
beta0fwhm = beta0
beta1fwhm = beta1
betamfwhm = betam
residfwhm = outresid
if options.do_pc1:
tp1fwhm = os.path.join(dirname,'tp1.mgh')
#cmd = 'mri_surf2surf --hemi '+hemi+' --s '+tpids[0]+' --sval '+measures[0]+' --src_type curv --tval '+tp1fwhm+' --trg_type mgh'
cmd = 'mri_convert '+measures[0]+' '+tp1fwhm
run_cmd(cmd,'mri_convert converting TP1 measure did not work?')
# create ouput (depending on type)
if options.resid > 0:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing residuals tp "+str(options.resid)+"\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_resid+fwhmstr+'.mgh')
# copy (potentially smoothed) residual map to base (final within-subject output)
cmd = 'cp '+residfwhm+' '+outname
run_cmd(cmd,'cp residuals failed?')
allout.append(outname)
if options.qcache is not None:
outtname= os.path.join(basesurfdir,hemi+'.'+options.name_resid+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
if options.do_pc1fit:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing PC1FIT\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_pc1fit+fwhmstr+'.mgh')
# compute symmetrized pc1fit change:
cmd = 'mris_calc -o '+outname+' '+beta1fwhm+' div '+beta0fwhm
run_cmd(cmd,'mris_calc compute pct. change (pc1fit) problem?')
cmd = 'mris_calc -o '+outname+' --label '+label+' '+outname+' mul 100'
run_cmd(cmd,'mris_calc compute pct. change (pc1fit) problem?')
allout.append(outname)
if options.qcache is not None:
outtname= os.path.join(basesurfdir,hemi+'.'+options.name_pc1fit+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
if options.do_spc:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing SPC\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_spc+fwhmstr+'.mgh')
# compute symmetrized pct change:
cmd = 'mris_calc -o '+outname+' '+beta1fwhm+' div '+betamfwhm
run_cmd(cmd,'mris_calc compute sym. pct. change (spc) problem?')
cmd = 'mris_calc -o '+outname+' --label '+label+' '+outname+' mul 100'
run_cmd(cmd,'mris_calc compute sym. pct. change (spc) problem?')
allout.append(outname)
if options.qcache is not None:
outtname= os.path.join(basesurfdir,hemi+'.'+options.name_spc+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
if options.do_pc1:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing PC1\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_pc1+fwhmstr+'.mgh')
# compute pct change:
cmd = 'mris_calc -o '+outname+' '+beta1fwhm+' div '+tp1fwhm
run_cmd(cmd,'mris_calc compute percent change (pc1) problem?')
cmd = 'mris_calc -o '+outname+' --label '+label+' '+outname+' mul 100'
run_cmd(cmd,'mris_calc compute percent change (pc1) problem?')
allout.append(outname)
if options.qcache is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_pc1+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
if options.do_rate:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing Rate\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_rate+fwhmstr+'.mgh')
# mul 1 to remove everything outside the label:
cmd = 'mris_calc -o '+outname+' --label '+label+' '+beta1fwhm+' mul 1'
run_cmd(cmd,'mris_calc masking with label (rate)')
allout.append(outname)
if options.qcache is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_rate+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
if options.do_avg:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" fwhm "+str(fwhm)+" computing AVG\n"
outname = os.path.join(basesurfdir,hemi+'.'+options.name_avg+fwhmstr+'.mgh')
# mul 1 to remove everything outside the label:
cmd = 'mris_calc -o '+outname+' --label '+label+' '+betamfwhm+' mul 1'
run_cmd(cmd,'mris_calc masking with label (avg)')
allout.append(outname)
if options.qcache is not None:
outtname = os.path.join(basesurfdir,hemi+'.'+options.name_avg+fwhmstr+'.'+options.qcache+'.mgh')
alltout.append(outtname)
# end iterate fwhm levels
# now we have all outputs at all fwhm levels
# if qcache, map them to target:
if options.qcache is not None:
# stack ALL measure files and map to template and unstack again:
if len(allout) > 0:
print "\n==============================================================================="
print "SUBJECT "+subjectid+" mapping all measure files to QCache (as stack)\n"
stack_map(subjectid,options.qcache,'mgh','mgh', allout , alltout ,dirname,hemi, options.jac)
# cleanup tmp dir:
shutil.rmtree(dirname)
if len(allout) > 0:
print 'You can look at the result with, e.g.:'
print ' tksurfer '+subjectid+' '+hemi+' pial -overlay '+allout[-1]
print
# end iterate subjects
# stuff on common qcache target:
if options.qcache is not None:
if options.isec_labels is not None:
print "\n==============================================================================="
print "Intersecting all lables in template QCache\n"
# intersect labels on common qcache target:
first=True
for label in tlabels:
if first:
cmd = 'cp '+label+' '+options.isec_labels
run_cmd(cmd,'cp first qcache label failed?')
first = False
else:
cmd = "mris_label_calc intersect "+options.isec_labels+' '+label+' '+options.isec_labels
run_cmd(cmd,'mris_label_calc intersect failed?')
# erode a little
tpial = os.path.join(subjectsdir,options.qcache,"surf",hemi+'.pial')
if not os.path.exists(tpial):
print 'ERROR: qtemplate pial '+str(tpial)+' does not exist!'
sys.exit(1)
cmd="mris_label_calc erode 2 "+options.isec_labels+" "+tpial+" "+options.isec_labels
run_cmd(cmd,'mris_label_calc erode failed?')
# stack maps on qcache geometry if wanted:
if options.stack_spc is not None or options.stack_pc1fit is not None or options.stack_pc1 is not None or options.stack_avg is not None or options.stack_rate is not None or options.stack_resid is not None:
print "\n==============================================================================="
print "Stacking all maps in template QCache\n"
fwhms = [ 0 ]
if options.fwhm is None:
fwhms = [ 0 , 5, 10, 15, 20, 25 ]
elif options.fwhm is not None:
fwhms = [ options.fwhm ]
for fwhm in fwhms:
fwhmstr = ''
if fwhm > 0 :
fwhmstr = '.fwhm'+str(fwhm)
# stack each type of measure for this fwhm:
if options.stack_resid is not None:
allin = " ".join(tresid[fwhm])
allstack=options.stack_resid
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache resid outputs at fwhm '+str(fwhm)+' failed?')
if options.stack_spc is not None:
allin = " ".join(tspc[fwhm])
allstack=options.stack_spc
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache spc outputs at fwhm '+str(fwhm)+' failed?')
if options.stack_pc1fit is not None:
allin = " ".join(tpc1fit[fwhm])
allstack=options.stack_pc1fit
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache pc1fit outputs at fwhm '+str(fwhm)+' failed?')
if options.stack_pc1 is not None:
allin = " ".join(tpc1[fwhm])
allstack=options.stack_pc1
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache pc1 outputs at fwhm '+str(fwhm)+' failed?')
if options.stack_avg is not None:
allin = " ".join(tavg[fwhm])
allstack=options.stack_avg
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache avg outputs at fwhm '+str(fwhm)+' failed?')
if options.stack_rate is not None:
allin = " ".join(trate[fwhm])
allstack=options.stack_rate
if allstack[-4:] == ".mgh": allstack = allstack[:-4]
allstack= allstack+fwhmstr+'.mgh'
cmd = 'mri_concat --i '+allin+' --o '+allstack
run_cmd(cmd,'mri_concat stacking qcache rate outputs at fwhm '+str(fwhm)+' failed?')
# always exit with 0 exit code
sys.exit(0)