#!/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)