#! /usr/bin/env perl

############################################################################
############################################################################
### Run CIVET processing pipeline within the PMP framework. 
###         
###   In the current version, CIVET requires a quarantine of the software
###   invoked and an environment file.
###                                       
###   Authors:                                  
###         Original pipeline scripts by: Jason Lerch 
###         Pre-modular versions of CIVET: Yasser Ad-Dab'bagh             
###         Modular CIVET: Oliver Lyttelton, J-Sebastian Muehlboeck, 
###                        Yasser Ad-Dab'bagh, Kelvin Mok, Claude Lepage
###         Project started in May 3, 2005       
###
###         Copyright Alan C. Evans
###         Professor of Neurology
###         McGill University
###
### For more information, have a look at our documentation at:
###    http://wiki.bic.mni.mcgill.ca/index.php/CIVET
###                                 
############################################################################
############################################################################

use strict;
use FindBin;
use Cwd qw( abs_path );
use Env qw( PATH );

# All modules that will be used in either case are declared here
use Getopt::Tabular;
use PMP::PMP;
use PMP::pbs;
use PMP::sge;
use PMP::spawn;
use PMP::Array;
use MNI::Startup;
use MNI::PathUtilities qw(split_path);
use MNI::FileUtilities qw(check_output_dirs check_output_path);
use MNI::DataDir;

use lib "$FindBin::Bin";
use MRI_Image;
use Processing_Pipeline_Main;

$PATH = "$FindBin::Bin/progs:${PATH}";

# Set interrupt handler (for cleaning of lock files)
$SIG{'INT'} = 'CLEANUP';
$SIG{'TERM'} = 'CLEANUP';

my $version = "1.1.12";
my $versionDate= "February, 2013";
my $usage = "\nUSAGE:\n$ProgramName -sourcedir <dir> -targetdir <dir> -prefix <prefix> [options] id1 id2 ... idn > <logfile> &\n 
ALTERNATIVE USAGE:\n
$ProgramName -sourcedir <dir> -targetdir <dir> -prefix <prefix> -id-file <idfile> [options]  > <logfile> &\n\n";

my $whatsnew = <<NEW;
    *************************************************************************
    New in $versionDate quarantine, version $version

    1 - Neck-cropping option on native image (-headheight)
    2 - Improved linear registration (now using normalized mutual information
        and with model_mask on last fitting stage only)
    3 - Improved non-linear registration (now with a headmask and with
        reduced weight parameter, 2mm steps and higher sampling)
    4 - Application of population model's head mask to reduce background
        noise in image for mincbet
    5 - Final pve classification is unmasked so it can recover any small
        bits of tissue missing from the brain mask
    6 - Improved masking of white matter mask using cortical_surface
        with a new model (on which the non-linear transform is applied)
    7 - Use smoothing from depth_potential instead of diffuse for 
        cortical thickness and surface registration (more accurate smoothing
        and no more need for BLAS and LAPACK). Note: 20mm with old code is 
        like 30mm with new code.
    8 - Improved filling of ventricles using an improved atlas (for
        white surface extraction)
    9 - Fixed WM/GM border by not allowing CSF/BG next to WM (under
        brainstem, mostly)
    10 - New ADNI population models
    11 - Generation of basic glim file as part of the CIVET QC
    12 - New verify figures for Laplacian fit and surface-surface
         intersections

    More to come soon:

    1 - improved medial cut through the corpus callosum.
    2 - bug fixed in N3 to be invariant to voxel size
    3 - high-resolution surfaces at 320k polygons.
    4 - classification and surface extraction in native space.
    5 - intermediate model for linear and non-linear volume registration.
    6 - Plug-in for alternate population templates
    7 - removal of self-intersections in .sm file for surface
        registration (Maxime Boucher)
    8 - Improved surface registration by geometric smoothing of
        the surface
    9 - Longitudinal surface registration and QC
    10 - Marching-cubes for initial white surface

    ********** From Nov-2010 quarantine, version 1.1.11 **********

    1 - option for user-defined surface parcellation atlas.
    2 - AAL surface parcellation atlas in models.
    3 - normalized cortical volumes maps on resampled surfaces
        (Lu Zhao and Maxime Boucher).

    ********** From Feb-2010 quarantine, version 1.1.10 **********

    1 - option for 3 Tesla.
    2 - options for reset-after and reset-to (PMP).
    3 - cortical area expansion/contraction maps on resampled surfaces.
    4 - new ICBM 152 non-linear model (40th generation, 2009),
        symmetric and asymmetric.
    5 - asymmetry maps for cortical thickness (now normalized)
    6 - asymmetry maps for position
    7 - basic lobar surface parcellation for quality control
        of surface registration
    8 - sum of absolute mean curvature
    9 - gyrification index for whole brain (in addition to by hemisphere)
    10 - GI and mean curvature on white, gray, mid surfaces

    ********** From Feb-14-2008 quarantine, version 1.1.9 **********

    1 - minc2-based quarantine.
    2 - new ICBM templates compatible with non-linear ICBM152 model.
    3 - asymmetry maps for cortical thickness.
    4 - cortical surfaces for hemispheres combined into single surface
        for convenience of visualization.
    5 - integrated VBM pipeline inside CIVET.
    6 - list of references produced by CIVET in targetdir.
    7 - really mask the background t1-image during classification to 
        avoid confusion between noisy background and csf.
    8 - fixed tlaplace cortical thickness for extents of volume in native space.
    9 - ANIMAL segmentation only available as a plug-in module.
    10 - resampled cortical surfaces in stereotaxic space.
    11 - total cerebral volume.

    ********** From Jul-09-2007 quarantine, version 1.1.8 **********

    1 - ANIMAL segmentation has been disabled as the new default.
    2 - recompute non-linear transform to linear symmetric icbm152 
        model when using the old ANIMAL parcellation.
    3 - simplification of PMP options and new -no-granular mode to
        submit all stages for one subject in one qsub job.
    4 - susceptibility artefact map now uses brain with cerebellum.
    5 - use a t1-native user-defined custom brain mask, if present,
        of the form SOURCE_DIR/PREFIX_ID_mask.mnc (or .gz) for
        linear registration and white matter masking.
    6 - use t1only for native brain mask (no need for multispectral
        mask in native).
    7 - run nu_correct in stereotaxic space on the interpolated 
        original native image, not on the interpolated nu_corrected
        native image.
    8 - init environment files are now in the quarantine base 
        directory (no longer in CIVET directory).

    ********** From previous Mar-30-2007 quarantine, version 1.1.7 **********

    1 - added mean curvature on native mid surfaces.
    2 - added gyrification index.
    3 - cortical thickness now always produced, with tlink method 
        and 20mm blurring kernel as defaults.
    4 - added percentages for tissue types in verify image.
    5 - much improved brain mask, remove -crop-neck option.

    ********** From previous Mar-20-2007 quarantine, version 1.1.6 **********

    1 - support for various interpolation methods from native to
        stereotaxic (linear, cubic, sinc; default is linear).
    2 - error metric for brain mask in native space (in verify image).
    3 - shuffle the PMP steps so that -resetAll actually works.
    4 - file name extension for surface maps is .sm, not .mnc.
        (use rename 's/\.mnc\$/\.sm/' */transforms/surfreg/*.mnc).
    5 - can use a 0mm blurring kernel for cortical thickness (no blurring).
    6 - improved cut through callosum to split hemispheres.
    7 - fixed bug in surface registration for infinite loop (on RedHat).
    8 - clear history in minc files to go around likely bug in netCDF 3.6.1.

    ********** From previous Feb-28-2007 quarantine, version 1.1.5 **********

    1 - surface registration and resampling of cortical thickness.
    2 - binary masks in byte format to save disk space.
    3 - specification of data-processing template for voxel size
        (0.50mm, 0.75mm, 1.0mm, 2.0mm, etc; default is 1.0mm).
    4 - check for irregular slice spacing in native images.
    5 - summary of options printed in logs directory of subject.
    6 - validation of input parameters for correct values.
    7 - possibility to easily rerun CIVET with different options
        for cortical thickness method and fwhm.
    8 - keep mid-surfaces in stereotaxic space
    9 - keep surface lobes parcellation at vertices

    ********** From previous Jan-24-2007 quarantine, version 1.1.4 **********

    1 - add surface segmentation based on ANIMAL\'s volume segmentation
        and calculation of cortex area for lobes (eventually will use
        surface registration with surface parcellation).
    2 - iterative correction to the mean and variance of the tissue types
        inside pve as the process converges (option -correct-pve)

    ********** From previous Jan-15-2007 quarantine, version 1.1.3 **********

    1 - add -nothreshold to mritoself for -spectral_mask. This affects only
        the multispectral mask, in very rare cases. If the native mask looks
        wrong and linear registration for no obvious reason, with -spectral_mask,
        then re-run your subject with this fix.
    2 - add -lsq6 to bestlinreg.pl and apply initial transform, if given, to
        mask, only if the mask is given too (does not apply to CIVET).

    ********** From previous Dec-20-2006 quarantine, version 1.1.2 **********

    1 - use the icbm152 mask for nu_correct instead of the avg305 mask,
        since subject is linearly registered to icbm152, not avg305.
    2 - transformation of surfaces to native space, with cortical thickness
        evaluated in native space.
    3 - quarantine based on minc-1.5 (upgraded versions of pcre, pcre++,
        CGAL,BOOST for compilation in 64 bits).
    4 - use icbm152 surface mask for verify image (red lines). Was avg305
        before.
    5 - use classified image from pve instead of classify_clean as input
        to cortical_mask for white matter masking.
    6 - improved mincbet for hyperintense voxels in t1 for masking meninges
    7 - verify.png now shows, in the first row, the brain mask that was 
        used for linear registration (linear registration can live with a
        mask that is slightly off, but not too wrong however)
    8 - make all.q the default queue instead of aces.q
    9 - new flag -lsq12 to CIVET for 12-param linear registration (default
        is still 9-param)
   10 - multispectral classification now produces the good result

    ********** From previous Nov-02-2006 quarantine, version 1.1.0 **********

    The major flaws in the Sep-12-2006 quarantine have been discovered and
    fixed and the former quarantine should no longer be used for any reason.
    The Sep-12-2006 quarantine would give too much white matter, thus bad 
    white surfaces and consequently wrong cortical thickness. All simulations 
    done with the Sep-12-2006 quarantine should be redone using the current 
    one.

    1 - nu_correct must be performed in Talairach space after linear
        registration, where a suitable brain mask is available.
    2 - multispectral brain mask now available in native space (was
        available only in Talairach space before)
    3 - fixed minc interface to mincbet to copy direction cosines of
        axes (was responsible for small rotation of brain in Talairach)
    4 - linear registration is more robust and is now really "best"
    5 - application of the non-linear transformation to the tag points
        used by classify_clean (essential for young children with small
        ventricles or AD subjects with large ventricles)

    There are several other minor enhancements that will not change the
    results but that will make your like easier:

    1 - verify.png image has an outline for the brain mask in Talairach 
        so that you can easily see if the linear registration has been
        successful
    2 - clean-up of the lock files when CIVET is interrupted
    3 - source and target directories may use relative paths
    4 - source files may or may not be gzipped. 

    *************************************************************************
NEW

my $help = <<HELP;

$ProgramName, version $version, released $versionDate.

    Takes any number of multi or single spectral input MINC volumes and
    extracts the cortical surfaces from them utilizing the PMP pipeline
    system. It then calculates cortical thickness at each vertex of the 
    produced cortical surfaces (non-linearly registered) using the t-link 
    metric (in both Talairach and native spaces). It can also produce ANIMAL 
    segmentations, symmetry analyses, regional thickness, surface areas 
    and volumes for brain lobes.

$whatsnew
HELP
Getopt::Tabular::SetHelp($help, $usage);

####################
# Argument handling:
####################

my $resetAll = undef;
my $resetFrom = undef;
my $resetAfter = undef;
my $resetTo = undef;
my $resetRunning = 1;

my $command = "printStatus";

# Default queuing defaults:

my %PMPconf = ( 'DEFAULT' => { 'type' => "spawn",
                               'maxqueued' => 10000,
                               'granularity' => 1,
                               'command' => undef,
                               'queue' => undef,
                               'hosts' => undef,
                               'opts' => undef },
                'MNIBIC'  => { 'type' => "sge",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "qsub",
                               'queue' => "all.q",
                               'hosts' => undef,
                               'opts' => undef },
                'ZLAB'    => { 'type' => "sge",
                               'maxqueued' => 100,
                               'granularity' => 1,
                               'command' => "qsub",
                               'queue' => "all.q",
                               'hosts' => undef,
                               'opts' => undef },
                'CLUMEQ'  => { 'type' => "pbs",     ## this was krylov
                               'maxqueued' => 100,
                               'granularity' => 1,
                               'command' => "msub",
                               'queue' => "brain",
                               'hosts' => undef,
                               'opts' => "-l ncpus=1" },
              'GUILLIMIN' => { 'type' => "pbs",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "msub",
                               'queue' => undef,
                               'hosts' => undef,
                               'opts' => "-l walltime=8:00:00 -l nodes=1:ppn=1" },
                'COLOSSE'  => { 'type' => "pbs",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "msub",
                               'queue' => "default",
                               'hosts' => undef,
                               'opts' => "-P eim-670-aa -l h_rt=6:00:00" },
                'RQCHP'   => { 'type' => "pbs",
                               'maxqueued' => 800,
                               'granularity' => 1,
                               'command' => "qsub",
                               'queue' => "qwork\@ms",
                               'hosts' => undef,
                               'opts' => "-l walltime=7:00:00" },
                'KISTI'   => { 'type' => "pbs",
                               'maxqueued' => 100,
                               'granularity' => 1,
                               'command' => "msub",   ## check
                               'queue' => undef,      ## check
                               'hosts' => undef,
                               'opts' => undef },
                'NIH'     => { 'type' => "pbs",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "qsub",
                               'queue' => "norm",
                               'hosts' => undef,
                               'opts' => "-l nodes=1:p2800" },
                'NUPPI'   => { 'type' => "sge",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "qsub",
                               'queue' => undef,    ## check
                               'hosts' => undef,
                               'opts' => undef },
                'JUDGE'   => { 'type' => "pbs",
                               'maxqueued' => 1000,
                               'granularity' => 1,
                               'command' => "msub",
                               'queue' => "common",
                               'hosts' => undef,
                               'opts' => "-l walltime=08:00:00" },
                'JUROPA'  => { 'type' => "pbs",      # should never run on Juropa
                               'maxqueued' => 1000,  # because ppn=16 not filled
                               'granularity' => 1,
                               'command' => "msub",
                               'queue' => "common",
                               'hosts' => undef,
                               'opts' => "-l walltime=08:00:00 -l nodes=1:ppn=1 -v tpt=2" } );

my $PMPtype = ( $ENV{'CIVET_JOB_SCHEDULER'} || "DEFAULT" );
my $PMPmaxQueued = $PMPconf{$PMPtype}{maxqueued};
my $PMPgranularity = $PMPconf{$PMPtype}{granularity};
my $PMPqueue = $PMPconf{$PMPtype}{queue};
my $PMPhosts = $PMPconf{$PMPtype}{hosts};
my $PMPopts = $PMPconf{$PMPtype}{opts};

############# Default models for pipeline

my $CIVETModel = "icbm152nl";
my $TemplateSize = "1.00";

### add avg64 model for cortical_surface, sphere models, 

my %CIVETmodels = ( 'icbm152nl_09a' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "mni_icbm152_t1_tal_nlin_asym_09a",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "mni_icbm152_t1_tal_nlin_asym_09a",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "mni_icbm152_t1_tal_nlin_asym_09a_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlasLeft'    => "$FindBin::Bin/models/surface_atlas.txt",
                      'SurfAtlasRight'   => "$FindBin::Bin/models/surface_atlas.txt",
                      'WhiteSurfMask'    => "$FindBin::Bin/models/Cerebellum_Ventricles_SubCortical_Mask.mnc",
                      'TagFileDir'       => MNI::DataDir::dir("classify"),
                      'TagFile'          => "ntags_1000_prob_90_nobg.tag",
                      'bgTagFile'        => "ntags_1000_bg.tag",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    },
                    'icbm152nl_09s' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "mni_icbm152_t1_tal_nlin_sym_09a",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "mni_icbm152_t1_tal_nlin_sym_09a",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "mni_icbm152_t1_tal_nlin_sym_09a_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlasLeft'    => "$FindBin::Bin/models/surface_atlas.txt",
                      'SurfAtlasRight'   => "$FindBin::Bin/models/surface_atlas.txt",
                      'WhiteSurfMask'    => "$FindBin::Bin/models/Cerebellum_Ventricles_SubCortical_Mask.mnc",
                      'TagFileDir'       => MNI::DataDir::dir("classify"),
                      'TagFile'          => "ntags_1000_prob_90_nobg.tag",
                      'bgTagFile'        => "ntags_1000_bg.tag",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    },
                    'icbm152nl' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "icbm_avg_152_t1_tal_nlin_symmetric_VI",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "icbm_avg_152_t1_tal_nlin_symmetric_VI",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "icbm_avg_152_t1_tal_nlin_symmetric_VI_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlasLeft'    => "$FindBin::Bin/models/surface_atlas.txt",
                      'SurfAtlasRight'   => "$FindBin::Bin/models/surface_atlas.txt",
                      'WhiteSurfMask'    => "$FindBin::Bin/models/Cerebellum_Ventricles_SubCortical_Mask.mnc",
                      'TagFileDir'       => MNI::DataDir::dir("classify"),
                      'TagFile'          => "ntags_1000_prob_90_nobg.tag",
                      'bgTagFile'        => "ntags_1000_bg.tag",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    },
                    'icbm152lin' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "icbm_avg_152_t1_tal_lin_symmetric",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "icbm_avg_152_t1_tal_lin_symmetric",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "icbm_avg_152_t1_tal_lin_symmetric_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlasLeft'    => "$FindBin::Bin/models/surface_atlas.txt",
                      'SurfAtlasRight'   => "$FindBin::Bin/models/surface_atlas.txt",
                      'WhiteSurfMask'    => "$FindBin::Bin/models/Cerebellum_Ventricles_SubCortical_Mask.mnc",
                      'TagFileDir'       => MNI::DataDir::dir("classify"),
                      'TagFile'          => "ntags_1000_prob_90_nobg.tag",
                      'bgTagFile'        => "ntags_1000_bg.tag",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    },
                    'ADNInl' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "mni_adni_t1w_tal_nlin_asym",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "mni_adni_t1w_tal_nlin_asym",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "mni_adni_t1w_tal_nlin_asym_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlasLeft'    => "$FindBin::Bin/models/surface_atlas.txt",
                      'SurfAtlasRight'   => "$FindBin::Bin/models/surface_atlas.txt",
                      'WhiteSurfMask'    => "$FindBin::Bin/models/Cerebellum_Ventricles_SubCortical_Mask-ADNI.mnc",
                      'TagFileDir'       => MNI::DataDir::dir("classify-ADNI"),
                      'TagFile'          => "mni_adni_t1w_tal_nlin_asym_ntags_1000_nobg.tag",
                      'bgTagFile'        => "mni_adni_t1w_tal_nlin_asym_ntags_1000_bg.tag",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    },
                  );

############# The status report filename:
my $statusReportFile = "CIVET_status_report.csv";

############# Options for how CIVET is to be run: the default will be to run
############# t1 data, cleaning tags, ANIMAL, the production of 82k 
############# polygon surfaces, no symmetry stages, and surface registration. 
my $inputType = "t1only";
my $correctPVE = 0;
my $maskType = "t1only";
my $interpMethod = "trilinear";
my $headheight = 0;
my $N3Distance = undef;
my $N3Damping = "1.0e-07";
my $lsqtype = "-lsq9";
my $Area_fwhm = "40";
my $Volume_fwhm = "40";
my $VBM = "noVBM";
my $VBM_fwhm = 8;
my $VBM_symmetry = "noSymmetry";
my $VBM_cerebellum = "Cerebellum";
my $surface = "SURFACE";
my $MaskBloodVessels = 0;
my @thickness = ("tlink","20");
my $ResampleSurfaces = 0;
my $MeanCurvature = 0;
my $CombineSurfaces = 0;
my @SurfaceAtlas = ();

my $animal = "noANIMAL";
my $AnimalAtlas = undef;
my $AnimalAtlasDir = undef;
my $AnimalNLRegDir = undef;
my $AnimalNLRegModel = undef;

############# Some directories that the user will need to specify
my $sourceDir = undef;
my $target = undef;
my $prefix = undef;
my $sourceSubDir = "noIdSubDir";
my $idTextFile = undef;

############# Options table
my @leftOverArgs;
my @argTbl = (
          ["Execution control", "section"],
          ["-spawn", "const", "DEFAULT", \$PMPtype, "Use the perl system interface to spawn jobs [default: use local host scheduler $PMPtype]"],
          ["-queue", "string", 1, \$PMPqueue, "Which queue to use", "<queue>"],
          ["-hosts", "string", 1, \$PMPhosts, "Colon separated list of hosts", "<hosts>"],
          ["-qopts", "string", 1, \$PMPopts, "Extra options to queuing system", "<opts>"],
          ["-no-granular|-granular", "boolean", 1, \$PMPgranularity, 
           "Granularity level for submission of jobs using queueing system."],
          ["-maxqueued", "string", 1, \$PMPmaxQueued, 
           "Maximum number of jobs that can be submitted at once.", "<val>"],

          ["File options", "section"],
          ["-sourcedir", "string", 1, \$sourceDir, "Directory containing the source files.", "<dir>"],
          ["-targetdir", "string", 1, \$target, "Directory where processed data will be placed.", "<dir>"],
          ["-prefix", "string", 1, \$prefix, "File prefix to be used in naming output files.", "<prefix>"],
          ["-id-subdir", "const", "IdSubDir", \$sourceSubDir, 
           "Indicate that the source directory contains sub-directories for each id"],
          ["-id-file", "string", 1, \$idTextFile, "A text file that contains all the subject id\'s (separated by space, tab, return or comma) that CIVET will run on.", "<file>"],

          ["Pipeline options", "section"],
          ["-template", "string", 1, \$TemplateSize, 
           "Define the template for image processing in stereotaxic space (0.50, 0.75, 1.00, 1.50, 2.00, 3.00, 4.00, 6.00).", "<val>"],
          ["-model", "string", 1, \$CIVETModel,
           "Define the model for image-processing: " .
           "\"icbm152nl_09s\" (MNI ICBM152 non-linear symmetric (2009)), " .
           "\"icbm152nl_09a\" (MNI ICBM152 non-linear asymmetric (2009)), " .
           "\"icbm152nl\" (MNI ICBM152 non-linear 6th generation), " .
           "\"icbm152lin\" (MNI ICBM152 linear), " .
           "\"ADNInl\" (MNI ADNI non-linear) ",
           "<model>"],
          ["-surface-atlas", "string", 2, \@SurfaceAtlas,
           "Define the atlas for surface parcellation ",
           "<atlas_left.txt> <atlas_right.txt>"],
          ["CIVET options", "section"],  
          ["-multispectral", "const", "multispectral", \$inputType, 
           "Use T1, T2 and PD native files for tissue classification."],
          ["-correct-pve|-no-correct-pve", "boolean", 1, \$correctPVE, 
           "Apply correction to the mean and variance of tissue types at pve iterations."],
          ["-spectral_mask", "const", "multispectral", \$maskType, 
           "Use T1, T2 and PD stereotaxic files for brain masking."],
          ["-interp", "string", 1, \$interpMethod,
           "Interpolation method from native to stereotaxic space (\"trilinear\", \"tricubic\", \"sinc\"",
           "<method>"],
          ["-headheight", "string", 1, \$headheight, 
           "head height in mm for neck cropping (suggested value: 170; default 0=none).", "<dist>"],
          ["-N3-distance", "string", 1, \$N3Distance, 
           "N3 spline distance in mm (suggested values: 200 for 1.5T scan; 50 for 3T scan).", "<dist>"],
          ["-N3-damping", "string", 1, \$N3Damping, 
           "N3 damping coefficient (lambda) (suggested values: 1.0e-07).", "<lambda>"],
          ["-lsq6", "const", "-lsq6", \$lsqtype,
           "use 6-parameter transformation for linear registration [default -lsq9]" ],
          ["-lsq12", "const", "-lsq12", \$lsqtype,
           "use 12-parameter transformation for linear registration [default -lsq9]" ],
          ["-no-surfaces", "const", "noSURFACE", \$surface, "don\'t build surfaces"],
          ["-mask-blood-vessels|-no-mask-blood-vessels", "boolean", 1, \$MaskBloodVessels,
           "mask blood vessels prior to white surface extraction"],
          ["-thickness", "string", 2, \@thickness,
           "compute cortical thickness and blur [tlink|tlaplace|tnear|tnormal] [kernel size in mm]"],
          ["-resample-surfaces|-no-resample-surfaces", "boolean", 1, \$ResampleSurfaces,
           "resample cortical surfaces"],
          ["-mean-curvature|-no-mean-curvature", "boolean", 1, \$MeanCurvature,
           "produce mean curvature maps on surfaces"],
          ["-area-fwhm", "string", 1, \$Area_fwhm, 
           "blurring kernel size in mm for resampled surface areas", "<fwhm>"],
          ["-volume-fwhm", "string", 1, \$Volume_fwhm, 
           "blurring kernel size in mm for resampled surface volumes", "<fwhm>"],
          ["-combine-surfaces|-no-combine-surfaces", "boolean", 1, \$CombineSurfaces,
           "combine left and right cortical surfaces"],

          ["VBM options", "section"],  
          ["-VBM", "const", "VBM", \$VBM, 
           "process VBM files for analysis [default -no-VBM]"],
          ["-no-VBM", "const", "noVBM", \$VBM, "don\'t process VBM files for analysis"],
          ["-VBM-fwhm", "string", 1, \$VBM_fwhm, 
           "blurring kernel size in mm for volume", "<fwhm>"],
          ["-VBM-symmetry", "const", "Symmetry", \$VBM_symmetry, 
           "run symmetry tools [default -no-VBM-symmetry]"],
          ["-no-VBM-symmetry", "const", "noSymmetry", \$VBM_symmetry, 
           "don\'t run symmetry tools"],
          ["-VBM-cerebellum", "const", "Cerebellum", \$VBM_cerebellum, 
           "keep cerebellum in VBM maps"],
          ["-no-VBM-cerebellum", "const", "noCerebellum", \$VBM_cerebellum, 
           "mask out cerebellum in VBM maps [default -VBM-cerebellum]"],

          ["ANIMAL options", "section"],  
          ["-animal", "const", "ANIMAL", \$animal, 
           "run volumetric ANIMAL segmentation and cortical surface lobe parcellation [default -no-animal]"],
          ["-no-animal", "const", "noANIMAL", \$animal, 
           "don\'t run volumetric ANIMAL segmentation and cortical surface lobe parcellation"],
          ["-symmetric_atlas", "const", "-symmetric_atlas", \$AnimalAtlas, 
           "Use linear symmetric atlas for ANIMAL segmentation (old way)"],
          ["-lobe_atlas", "const", "-lobe_atlas", \$AnimalAtlas, 
           "Use non-linear lobe atlas for ANIMAL segmentation (new way)"],
          ["-animal-atlas-dir", "string", 1, \$AnimalAtlasDir, 
           "Directory containing the segmentation atlas", "<dir>"],

          ["Pipeline control", "section"],
          ["-run", "const", "run", \$command, "Run the pipeline."],
          ["-status-from-files", "const", "statusFromFiles", \$command, "Compute pipeline status from files"],
          ["-print-stages", "const", "printStages", \$command, "Print the pipeline stages."],
          ["-print-status", "const", "printStatus", \$command, "Print the status of each pipeline."],
          ["-make-graph", "const", "makeGraph", \$command, "Create dot graph file."],
          ["-make-filename-graph", "const", "makeFilenameGraph", \$command, "Create dot graph of filenames."],
          ["-print-status-report", "const", "printStatusReport", \$command, "Writes a CSV status report to file in cwd."],

          ["Stage Control", "section"],
          ["-reset-all", "const", "resetAll", \$resetAll, "Start the pipeline from the beginning."],
          ["-reset-from", "string", 1, \$resetFrom, "Restart from the specified stage.", "<stage_name>"],
          ["-reset-after", "string", 1, \$resetAfter, "Restart after the specified stage.", "<stage_name>"],
          ["-reset-to", "string", 1, \$resetTo, "Run up to and including the specified stage.", "<stage_name>"],
          ["-reset-running|-no-reset-running", "boolean", 1, \$resetRunning, "Restart currently running jobs."],
          );

GetOptions(\@argTbl, \@ARGV, \@leftOverArgs) or die "\n";


############# Basic usage
my @dsids;

if ($idTextFile) {
    open (IDTEXTFILE, "$idTextFile") or die ("Cannot open '$idTextFile': $!");
    # read the whole text file into one string
    my $idstext = "";
    while (my $idline = <IDTEXTFILE>) {
    $idstext .= $idline;
    }
    close (IDTEXTFILE) or die ("Cannot close '$idTextFile': $!");
    # split the string on whitespace (\s) or comma
    @dsids = split(/[\s,]+/, $idstext);
} 
else {
    @dsids = @leftOverArgs or die $usage;
}

unless ($prefix && $target && $sourceDir) {
    die "\n\n*******ERROR********: \n     You must specify -prefix, -targetdir, and -sourcedir \n********************\n\n\n";
}
$target =~ s#/+$##;      # remove trailing / at end of directory name, if any
$target = abs_path( $target );
$sourceDir = abs_path( $sourceDir );

############# Override default PMP options based on command line options
$PMPconf{$PMPtype}{maxqueued} = $PMPmaxQueued;
$PMPconf{$PMPtype}{granularity} = $PMPgranularity;
$PMPconf{$PMPtype}{queue} = $PMPqueue;
$PMPconf{$PMPtype}{hosts} = $PMPhosts;
$PMPconf{$PMPtype}{opts} = $PMPopts;

############# Set no file buffering for stdout (buffer is printed every 1 line)
$| = 1;

############# Print the CIVET options list and related error messages
my $DATE = `date`;
chomp( $DATE );
my $UNAME = `uname -s -n -r`;
chomp( $UNAME );
print "\n\n* Pipeline started at $DATE on $UNAME \n";
print "\n* This is CIVET $version, $versionDate \n";
print "\n* CIVET Command line is:\n       $0 @ARGV \n";
print "\n* The source directory is: '$sourceDir' \n";
print "* The target directory is: '$target' \n";
print "* The prefix is: '$prefix' \n";
print "* The PMP class is: '$PMPconf{$PMPtype}{type}' \n";
print "* The $PMPtype queue type is: '$PMPconf{$PMPtype}{queue}' \n";
print "* The $PMPtype batch host(s) is/are: '$PMPconf{$PMPtype}{hosts}' \n";
print "* The $PMPtype global options are: '$PMPconf{$PMPtype}{opts}' \n";

############# The -like template (0.50, 0.75, 1.00, 1.50, 2.00, 3.00, 4.00, 6.00)

my $Template = "$CIVETmodels{$CIVETModel}{TemplateDir}/" .
               "$CIVETmodels{$CIVETModel}{TemplateModel}_${TemplateSize}mm.mnc";
print "* Template for image-processing is:\n  $Template \n";

############# Define registration targets and directories and print them

print "* Using $CIVETModel model for data-processing pipeline:\n";

my $regLinModelDir = $CIVETmodels{$CIVETModel}{RegLinDir};
my $regNLModelDir = $CIVETmodels{$CIVETModel}{RegNLDir};

print "* Registration linear model directory is:    \n  $regLinModelDir \n";
print "* Registration non-linear model directory is:\n  $regNLModelDir \n";

my $regLinModel = $CIVETmodels{$CIVETModel}{RegLinModel};
my $regNLModel = $CIVETmodels{$CIVETModel}{RegNLModel};

print "* Registration linear model is:    \n   $regLinModel \n";
print "* Registration non-linear model is:\n   $regNLModel \n";

############# Define intermediate registration targets and directories and print them
my $intermediate_model = undef;

############# Define Surface registration targets and directories and print them

my $SurfRegModelDir = $CIVETmodels{$CIVETModel}{SurfRegModelDir};
my $SurfRegModel = $CIVETmodels{$CIVETModel}{SurfRegModel};
my $SurfRegDataTerm = $CIVETmodels{$CIVETModel}{SurfRegDataTerm};

print "* Surface registration model directory is: \n  ${SurfRegModelDir} \n";
print "* Surface registration model is:    \n   ${SurfRegModel} \n";
print "* Surface registration dataterm is: \n   ${SurfRegDataTerm} \n";

############# Choose surface parcellation atlas if surfaces are to be resampled.

if( defined @SurfaceAtlas ) {
  if( !$ResampleSurfaces ) {
    die "You specified a surface parcellation atlas but did not ask to " .
        "resample the surfaces. Use -resample-surfaces option.\n";
  }
  if( defined $SurfaceAtlas[0] && defined $SurfaceAtlas[1] ) {
    if( -e $SurfaceAtlas[0] && -e $SurfaceAtlas[1] ) {
      $CIVETmodels{$CIVETModel}{SurfAtlasLeft} = $SurfaceAtlas[0];
      $CIVETmodels{$CIVETModel}{SurfAtlasRight} = $SurfaceAtlas[1];
    } else {
      if( -e "$CIVETmodels{$CIVETModel}{SurfRegModelDir}/$SurfaceAtlas[0]" && 
          -e "$CIVETmodels{$CIVETModel}{SurfRegModelDir}/$SurfaceAtlas[1]" ) {
        $CIVETmodels{$CIVETModel}{SurfAtlasLeft} = "$CIVETmodels{$CIVETModel}{SurfRegModelDir}/$SurfaceAtlas[0]";
        $CIVETmodels{$CIVETModel}{SurfAtlasRight} = "$CIVETmodels{$CIVETModel}{SurfRegModelDir}/$SurfaceAtlas[1]";
      } else {
        die "You specified a surface parcellation atlas but the atlas " .
            "cannot be found. Please check filename for atlas.\n";
      }
    }
  }
}

print "* Surface parcellation model (left) is: \n  $CIVETmodels{$CIVETModel}{SurfAtlasLeft}\n";
print "* Surface parcellation model (right) is: \n  $CIVETmodels{$CIVETModel}{SurfAtlasRight}\n";

############# Choose ANIMAL model.

if( $animal eq "ANIMAL" ) {
  if( ( defined $AnimalAtlas ) && 
      ( $AnimalAtlas eq "-symmetric_atlas" || $AnimalAtlas eq "-lobe_atlas" ) ) {
    if( !( defined $AnimalAtlasDir ) ) {
      # Try to get the default one in the quarantine, if it has been installed.
      $AnimalAtlasDir = MNI::DataDir::dir("ANIMAL_INSECT");
    }
    if( !( -e $AnimalAtlasDir ) ) {
      print "Warning: Required ANIMAL atlas directory $AnimalAtlasDir does not exist.\n";
      $animal = "noANIMAL";
    } else {
      if( $AnimalAtlas eq "-symmetric_atlas" ) {
        $AnimalAtlasDir .= "/icbm152-symmetric-v1.0/";
        $AnimalNLRegModel = "icbm_avg_152_t1_tal_lin_symmetric";
        $AnimalNLRegDir = MNI::DataDir::dir("mni-models");
      } elsif( $AnimalAtlas eq "-lobe_atlas" ) {
        $AnimalAtlasDir .= "/icbm152-lobes-v1.1/";
        $AnimalNLRegModel = "icbm_avg_152_t1_tal_nlin_symmetric_VI";
        $AnimalNLRegDir = MNI::DataDir::dir("mni-models");
      }
      if( !( -e $AnimalAtlasDir ) ) {
        print "Warning: Required ANIMAL atlas directory $AnimalAtlasDir does not exist.\n";
        $animal = "noANIMAL";
      }
    }
  } else {
    print "Warning: You must specify an atlas to use ANIMAL segmentation.\n";
    $animal = "noANIMAL";
  }
}

if( $animal eq "ANIMAL" ) {
  $CIVETmodels{$CIVETModel}{AnimalAtlas} = $AnimalAtlas;
  $CIVETmodels{$CIVETModel}{AnimalAtlasDir} = $AnimalAtlasDir;
  $CIVETmodels{$CIVETModel}{AnimalNLRegDir} = $AnimalNLRegDir;
  $CIVETmodels{$CIVETModel}{AnimalNLRegModel} = $AnimalNLRegModel;
  print "* ANIMAL atlas is $AnimalAtlas \n";
  print "* ANIMAL atlas directory is: \n   ${AnimalAtlasDir} \n";
  print "* ANIMAL non-linear registration model directory is: \n  ${AnimalNLRegDir} \n";
  print "* ANIMAL non-linear registration model is:    \n   ${AnimalNLRegModel} \n";
  print "\n";
  print "*             -WARNING- -WARNING- -WARNING-\n\n";
  print "* We recommend that you no longer use the ANIMAL segmentation in CIVET\n";
  print "* until the new lobe segmentation based on the iterative non-linear \n";
  print "* ICBM152 model is fully validated and released. Current results based\n";
  print "* on the linear symmetric model are inaccurate and should not be used.\n";
}


############# Print the Pipeline and Stage Control commands

if( !($command eq "run") ) {
	$resetRunning = 0;
}
print "\n\n* Pipeline Control command is: '$command'";

print "\n* Stage control commands are: " .
      ( $resetAll ) ? ( "reset all" ) : (
        ( $resetFrom ) ? ( "reset from '$resetFrom'" ) : (
          ( $resetAfter ) ? ( "reset after '$resetAfter'" ) : (
            ( $resetTo ) ? ( "reset to $resetTo" ) : ( "reset all" )
          )
        ) 
      ) . " and " .
      ( $resetRunning ) ? ( "'reset-running'" ) : ( "'no-reset-running'" );


print "\n\n\n* Data-set Subject ID(s) is/are: '@dsids'\n\n\n";

system("mkdir -p ${target}") if (! -d ${target});


############# An array to store the pipeline definitions for each subject
my $pipes = PMP::Array->new();

#########################################################
# Set up the pipeline output directories for each subject
#########################################################
foreach my $dsid (@dsids) {

    ##### Create image object. #####

    # depending on the two options your files can be in on target
    # or there can be subdirs for every subject
    my $Source_Base;
    if ($sourceSubDir eq "noIdSubDir") {
      $Source_Base = "${sourceDir}/";
    } else {
      $Source_Base = "${sourceDir}/${dsid}/";
    }

    my $image = MRI_Image->new( $version,
                                $Source_Base,
                                $target,
                                $prefix,
                                $dsid,
                                $inputType,
                                $correctPVE,
                                $maskType,
                                $interpMethod,
                                $headheight,
                                $MaskBloodVessels,
                                $N3Distance,
                                $N3Damping,
                                $lsqtype,
                                $surface,
                                \@thickness,
                                $ResampleSurfaces,
                                $MeanCurvature,
                                $Area_fwhm,
                                $Volume_fwhm,
                                $CombineSurfaces,
                                $VBM,
                                $VBM_fwhm,
                                $VBM_symmetry,
                                $VBM_cerebellum,
                                $animal,
                                $Template,
                                \$CIVETmodels{$CIVETModel} );

    ##### Create references once, for first subject only #####

    if( $dsid eq $dsids[0] ) {
      $image->make_references( $target );
    }

    ##### Create the pipeline environment. #####

    my $pipeline;

    if ($PMPconf{$PMPtype}{type} eq "spawn") {
        $pipeline = PMP::spawn->new();
        $PMPgranularity = 0;
    } elsif ($PMPconf{$PMPtype}{type} eq "sge") {
        $pipeline = PMP::sge->new();
        $pipeline->setCommand($PMPconf{$PMPtype}{command});
        $pipeline->setQueue($PMPconf{$PMPtype}{queue});
        $pipeline->setHosts($PMPconf{$PMPtype}{hosts}) if( defined $PMPconf{$PMPtype}{hosts} );
        $pipeline->setQueueOptions($PMPconf{$PMPtype}{opts}) if( defined $PMPconf{$PMPtype}{opts} );
	$pipeline->setPriorityScheme("later-stages");
    } elsif ($PMPconf{$PMPtype}{type} eq "pbs") {
        $pipeline = PMP::pbs->new();
        $pipeline->setCommand($PMPconf{$PMPtype}{command});
        $pipeline->setQueue($PMPconf{$PMPtype}{queue});
        $pipeline->setHosts($PMPconf{$PMPtype}{hosts}) if( defined $PMPconf{$PMPtype}{hosts} );
        $pipeline->setQueueOptions($PMPconf{$PMPtype}{opts}) if( defined $PMPconf{$PMPtype}{opts} );
        $pipeline->setPriorityScheme("later-stages");
    }

    # set some generic pipeline options
    $pipeline->name($dsid);
    $pipeline->debug(0);
    $pipeline->statusDir("${target}/${dsid}/logs");
    
##################################### ################# #################
# This used to be the definition of stages... Now it is in modules       
##################################### ################# #################

    CIVET_Main::create_pipeline(
                        \$pipeline,
                        \$image,
                        \$CIVETmodels{$CIVETModel},
                        $intermediate_model, 
                        $Template
                        );    

####################################################################################
# Pipeline management
####################################################################################

############# Set the status of the stages.

    if ($resetAll && ($resetAll eq "resetAll") ) {
        # Restart all stages
        $pipeline->resetAll();
    } else {
        # restart from a given stage;
        $pipeline->resetFromStage($resetFrom) if( $resetFrom );
        $pipeline->resetAfterStage($resetAfter) if( $resetAfter );
        $pipeline->subsetToStage($resetTo) if( $resetTo );
      
        # Rerun any failures from a previous run of this subjects pipe
        $pipeline->resetFailures();
        # Reset running jobs if that's what the user wants.
        $pipeline->resetRunning() if ($resetRunning);
    }

    $pipeline->updateStatus();
 
############# Add this pipe to our happy array of pipes
    
    $pipes->addPipe($pipeline);
    
} ####################################### end of foreach


############# Now run whatever it is that the user wanted done

if ($command eq "printStatus" ) {
    $pipes->printUnfinished();
}

elsif ($command eq "printStatusReport") {
    $pipes->printStatusReport($statusReportFile);
}

elsif ($command eq "statusFromFiles") {
    $pipes->updateFromFiles();
    $pipes->printUnfinished();
}

elsif ($command eq "printStages") {
    $pipes->printStages();
}

elsif ($command eq "makeGraph") {
    $pipes->createDotGraph("dependency-graph.dot");
}

elsif ($command eq "makeFilenameGraph") {
    $pipes->createFilenameDotGraph("filename-dependency-graph.dot","${target}/$dsids[0]/");
}

elsif ($command eq "run") { 
    # register all the programs
    $pipes->registerPrograms();
    $pipes->maxQueued($PMPmaxQueued);
    $pipes->setGranularity($PMPgranularity);
    $pipes->run();
}

else {
    print "huh? Grunkle little gnu, grunkle\n";
}

############# Voila!! #############

# Set interrupt handler (for cleaning of lock files)

sub CLEANUP {
  print "\nCaught signal - cleaning up all lock files\n\n";
  $pipes->cleanup();
  exit(1);
}