#! /bin/csh -f

#
# isxconcat-sess
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: greve $
#    $Date: 2012/11/27 21:03:24 $
#    $Revision: 1.60.2.1 $
#
# Copyright (C) 2002-2007,
# The General Hospital Corporation (Boston, MA). 
# All rights reserved.
#
# Distribution, usage and copying of this software is covered under the
# terms found in the License Agreement file named 'COPYING' found in the
# FreeSurfer source code root directory, and duplicated here:
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferOpenSourceLicense
#
# General inquiries: freesurfer@nmr.mgh.harvard.edu
# Bug reports: analysis-bugs@nmr.mgh.harvard.edu
#


set VERSION = '$Id: isxconcat-sess,v 1.60.2.1 2012/11/27 21:03:24 greve Exp $'

set inputargs = ($argv);

set outtop = ();
set analysis = ();
set flac = ();
set fsd  = ();
set contrastlist = ();

set conmaplist = ();
set anamaplist = (meanfunc raw.fsnr fsnr rho1mn)
set RegDOF = 6;
set regfile = ()
set interp = trilin
set fwhm = ();
set paireddiff = 0;
set pairedavg  = 0;
set mean = 0;
set ToVol = 0;
set ToSurf = 0;
set DoPercent = 0;
set DoAllContrasts = 0;
set DoOSGM = 0;
set DoWRFX = 0;
set DoRFX = 0;
set DoFFX = 0;
set SplitMVC = 1; # Split multivariate contrasts
set DoResample = 1;
set GetHRF = ();
set DoAnats = 0;

# for volume-based only
set talres = 2;
set talxfm = ();
set DoMorph = 0;

# for surface-based only
set hemilist = ();
set projfrac = (); # The default will be 0.5, set below
set DoProjFracAvg = 0;
set trgsubject = fsaverage; 
set fsgd = ();

set sessargs = ();
set tmpdir = ();
set cleanup = 1;
set PrintHelp = 0;

if($#argv == 0)  goto usage_exit;
set n = `echo $argv | grep -e version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif
set n = `echo $argv | grep -e -help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;  
  goto usage_exit;
endif

set SessList = `getsesspath $inputargs`;
if($status || $#SessList == 0) then
  echo "ERROR: finding sessions"
  echo "       $SessList" 
  exit 1;
endif

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:


# If the output exists, check consistency with group
set grouplistfile = $outtop/grouplist.txt 
if( -e $grouplistfile) then
  set grouplist = (`cat $grouplistfile`)
  if($#grouplist != $#SessList) then
    echo "ERROR: there is an inconsistency between the sessions already"
    echo "processed in $outtop and the current list of sesssions."
    echo "The previous group had $#grouplist sessions."
    echo "The current group has $#SessList sessions."
    echo "There are three things that you can do:"
    echo "  1. Specify a different output directory for the new group"
    echo "  2. Re-run with consistent group definition"
    echo "  3. Delete the output directory and re-run"
    exit 1;
  endif
  @ nth = 0
  foreach sess ($SessList)
    @ nth = $nth + 1
    set newsessid = `basename $sess`;
    set oldsessid = $grouplist[$nth];
    if($newsessid != $oldsessid) then
      echo "ERROR: there is an inconsistency between the sessions already"
      echo "processed in $outtop and the current list of sessions."
      echo "Session number $nth is $oldsessid in the previous analysis and "
      echo "  $newsessid in the current analysis"
      echo "There are three things that you can do:"
      echo "  1. Specify a different output directory for the new group"
      echo "  2. Re-run with consistent group definition"
      echo "  3. Delete the output directory and re-run"
      exit 1;
    endif
  end
endif

set outdir = $outtop/$analysis
mkdir -p $outdir/log

##### Create a log file ######
set LF = $outdir/log/isxconcat.log
if(-e $LF) mv $LF $LF.old
echo $LF

echo "isxconcat-sess log file" >> $LF
echo $VERSION >> $LF
date          >> $LF
uname -a      >> $LF
pwd           >> $LF
echo $0       >> $LF
echo $inputargs  >> $LF
id            >> $LF
echo GetHRF = $GetHRF >> $LF

set StartTime = `date`;

if($#flac == 0) cp $analysis/analysis.* $outdir

# Go thru each session, get DOF, etc
set sessinfofile = $outtop/sess.info.txt
set subjlistfile = $outtop/subjectlist.txt
set grouplistfiletmp = $outtop/grouplist.tmp.txt 
rm -f $sessinfofile $subjlistfile $grouplistfiletmp
@ DOF = 0;
foreach sess ($SessList)
  set sessid = `basename $sess`;
  echo $sessid >> $grouplistfiletmp

  set doffile = $sess/$fsd/$analysis/dof
  if(! -e $doffile) then
    echo "ERROR: cannot find $doffile" 
    exit 1;
  endif
  @ dofsess = `cat $doffile`;
  @ DOF = $DOF + $dofsess;

  set bstampfile = $sess/$fsd/$analysis/fs.build-stamp.txt
  if(-e $bstampfile) then
    set bstamp = `cat $bstampfile`
    set mdate = `stat -c %y $bstampfile | awk '{print $1}'`
  else
    set bstamp = ();
    set mdate = ();
  endif

  set srcsubject = `cat $sess/subjectname`
  echo "$sessid $srcsubject $dofsess $mdate $bstamp" >> $sessinfofile
  echo "$srcsubject" >> $subjlistfile
end
set ffxdofdat = $outdir/ffxdof.dat
echo $DOF > $ffxdofdat
set sess0 = $sess[1];
if(! -e $grouplistfile) cp $grouplistfiletmp $grouplistfile
rm -f $grouplistfiletmp

#-------------------------------------------------
foreach hemi ($hemilist)

  if($IsNative) then
    set SpStr  = $hemi.
    set SpStr2 = .$hemi
  else
    set SpStr = ()
    set SpStr2 = ()
  endif

  # Concat the masks
  set masksvol = $outdir/"$SpStr"masks.$FSF_OUTPUT_FORMAT
  set cmd = (isxconcat0-sess $sessargs -m mask -o $masksvol)
  if($#flac == 0) set cmd = ($cmd -a $analysis)
  if($#flac)      set cmd = ($cmd -flac $flac -fsd $fsd)
  if($#tmpdir)    set cmd = ($cmd -tmp $tmpdir);
  if($#fsgd) set cmd = ($cmd -fsgd $fsgd) # only need to do this once
  if($DoResample) then
    set cmd = ($cmd -interp nearest  -reg $regfile)
    if($hemi == mni305) then
      if(! $DoMorph) set cmd = ($cmd -talres $talres)
      if($DoMorph)   set cmd = ($cmd -morph)
      if($#talxfm)   set cmd = ($cmd -talxfm $talxfm)
    else
      set cmd = ($cmd -subject $trgsubject)
      set cmd = ($cmd -hemi $hemi)
      if($#projfrac)     set cmd = ($cmd -projfrac $projfrac)
      if($DoProjFracAvg) set cmd = ($cmd -projfrac-avg)
    endif
  endif
  echo "cd `pwd`" |& tee -a $LF
  echo $cmd       |& tee -a $LF
  $cmd            |& tee -a $LF
  if($status) exit 1;

  # Make a single mask = 0 unless all input masks are != 0
  set maskvol = $outdir/"$SpStr"mask.$FSF_OUTPUT_FORMAT
  set cmd = (mri_binarize --i $masksvol --frame-and --o $maskvol)
  echo $cmd       |& tee -a $LF
  $cmd            |& tee -a $LF
  if($status) exit 1;

  if($hemi == mni305) then
    # Take intersection of single mask and subcort mask
    set maskvol2 = $outdir/"$SpStr"subcort.mask.$FSF_OUTPUT_FORMAT
    set cmd = (mri_mask $maskvol $subcortmask $maskvol2)
    echo $cmd       |& tee -a $LF
    $cmd            |& tee -a $LF
    if($status) exit 1;
  endif

  # Go thru the maps in the contrasts
  foreach contrast ($contrastlist)
    set outcondir = $outdir/$contrast
    mkdir -p $outcondir

    # copy the contrast.mat file for reference
    if($#flac == 0) then
      # Might not exist for some reason
      set cmat = $analysis/$contrast.mat
      if(-e $cmat)  cp $cmat $outcondir/$contrast.lowerlevel.mat
    endif

    # Get number of variates
    set ces = `stem2fname $sess0/$fsd/$analysis/$contrast/ces`
    set tmpfile = /tmp/tmp.isxconcat-sess.$$.nframes
    mri_info --o $tmpfile --nframes $ces
    set J = (`cat $tmpfile`);
    rm $tmpfile
    echo $contrast J = $J |& tee -a $LF

    foreach mapname ($conmaplist)

      set outvol = $outcondir/"$SpStr"$mapname.$FSF_OUTPUT_FORMAT
      set cmd = (isxconcat0-sess -m $mapname -o $outvol $sessargs -c $contrast)
      if($#flac == 0) set cmd = ($cmd -a $analysis)
      if($#flac)      set cmd = ($cmd -flac $flac -fsd $fsd)
      if($SplitMVC) set cmd = ($cmd -split $tmpdir);
      if($#tmpdir)  set cmd = ($cmd -tmp $tmpdir);
      if($paireddiff) then
        if($mapname == cesvar || $mapname == cesvarpct) then
          set cmd = ($cmd -paired-sum) # Prob not right for var
        else
          set cmd = ($cmd -paired-diff)
        endif
      endif
      if($pairedavg)  set cmd = ($cmd -paired-avg)
      if($DoResample) then
        set cmd = ($cmd -interp $interp -reg $regfile)
        if($#fwhm)    set cmd = ($cmd -fwhm $fwhm);
        if($hemi == mni305) then
          if(! $DoMorph) set cmd = ($cmd -talres $talres)
          if($DoMorph)   set cmd = ($cmd --morph)
          if($#talxfm)   set cmd = ($cmd -talxfm $talxfm)
        else
          set cmd = ($cmd -hemi $hemi)
          set cmd = ($cmd -subject $trgsubject)
          if($#projfrac)      set cmd = ($cmd -projfrac $projfrac)
          if($DoProjFracAvg) set cmd = ($cmd -projfrac-avg)
        endif
      endif
      echo "cd `pwd`" |& tee -a $LF
      echo $cmd       |& tee -a $LF
      $cmd            |& tee -a $LF
      if($status) exit 1;

     end # Loop over maps

     if($DoOSGM) then
       set zp = ()
       if($J > 1) set zp = .000
       if(! $DoPercent ) then
         set y       = $outcondir/"$SpStr"ces$zp.$FSF_OUTPUT_FORMAT
         set yffxvar = $outcondir/"$SpStr"cesvar$zp.$FSF_OUTPUT_FORMAT
       else
         set y       = $outcondir/"$SpStr"cespct.$FSF_OUTPUT_FORMAT
         set yffxvar = $outcondir/"$SpStr"cesvarpct.$FSF_OUTPUT_FORMAT
       endif

       if($DoWRFX) then
         if(! $DoPercent) set glmdir = $outcondir/osgm.wrfx"$SpStr2"
         if($DoPercent) set glmdir = $outcondir/osgm.wrfx.pct"$SpStr2"
         set cmd = (mri_glmfit --osgm --y $y --glmdir $glmdir --nii)
         set cmd = ($cmd --mask $maskvol --wls $yffxvar --kurtosis)
         if($ToSurf) set cmd = ($cmd --surf $trgsubject $hemi --cortex)
         echo $cmd |& tee -a $LF
         $cmd |& tee -a $LF
         if($status) exit 1;
       endif

       if($DoRFX) then
         if(! $DoPercent) set glmdir = $outcondir/osgm.rfx"$SpStr2"
         if($DoPercent) set glmdir = $outcondir/osgm.rfx.pct"$SpStr2"
         set cmd = (mri_glmfit --osgm --y $y --glmdir $glmdir --kurtosis)
         set cmd = ($cmd --mask $maskvol --nii --kurtosis)
         if($ToSurf) set cmd = ($cmd --surf $trgsubject $hemi --cortex)
         echo $cmd |& tee -a $LF
         $cmd |& tee -a $LF
         if($status) exit 1;
       endif

       if($DoFFX) then
         if(! $DoPercent) set glmdir = $outcondir/osgm.ffx"$SpStr2"
         if($DoPercent) set glmdir = $outcondir/osgm.ffx.pct"$SpStr2"
         set cmd = (mri_glmfit --osgm --y $y --glmdir $glmdir --kurtosis)
         set cmd = ($cmd --mask $maskvol --yffxvar $yffxvar)
         set cmd = ($cmd --ffxdofdat $ffxdofdat --nii)
         if($ToSurf) set cmd = ($cmd --surf $trgsubject $hemi --cortex)
         echo $cmd |& tee -a $LF
         $cmd |& tee -a $LF
         if($status) exit 1;
        endif

      endif # OSGM

   end # Loop over contrasts

  # Analysis-level maps
  foreach mapname ($anamaplist)
    set outvol = $outdir/"$SpStr"$mapname.$FSF_OUTPUT_FORMAT
    set cmd = (isxconcat0-sess -m $mapname -o $outvol $sessargs);
    if($#flac == 0) set cmd = ($cmd -a $analysis)
    if($#flac)      set cmd = ($cmd -flac $flac -fsd $fsd)
    if($#tmpdir) set cmd = ($cmd -tmp $tmpdir);
    # What about paired diff for these?
    if($paireddiff && $mapname != h) then
      if($mapname == rvar) then
        set cmd = ($cmd -paired-sum)
      else
        set cmd = ($cmd -paired-diff)
      endif
    endif
    if($DoResample) then
      set cmd = ($cmd -interp $interp -reg $regfile )
      if($#fwhm) set cmd = ($cmd -fwhm $fwhm);
      if($hemi == mni305) then
        if(! $DoMorph) set cmd = ($cmd -talres $talres)
        if($DoMorph)   set cmd = ($cmd --morph)
        if($#talxfm)   set cmd = ($cmd -talxfm $talxfm)
      else
        set cmd = ($cmd -hemi $hemi)
        set cmd = ($cmd -subject $trgsubject)
        if($#projfrac)      set cmd = ($cmd -projfrac $projfrac)
        if($DoProjFracAvg) set cmd = ($cmd -projfrac-avg)
      endif
    endif
    if($mapname == h) then
      set cmd = ($cmd -gmean $#SessList)
      echo ""
      echo "Constructing the group HRF. This may take a while. If "
      echo "you do not want the HRF, re-run with -no-hrf"
      echo ""
      # It would be better to set the XtX to something like the identity
      if($ToVol)  cp $SessList[1]/$fsd/$analysis/h.dat $outdir/"$SpStr"h.dat
      if($ToSurf) cp $SessList[1]/$fsd/$analysis/h.dat $outdir/"$SpStr"h.dat
    endif
    echo "cd `pwd`" |& tee -a $LF
    echo $cmd       |& tee -a $LF
    $cmd            |& tee -a $LF
    if($status) exit 1;
  end # mapname

end # Loop over hemi


if($DoAnats) then
  foreach hemi ($hemilist)
    if($hemi != mni305) then
      foreach mapname (thickness curv sulc)
        set outvol = $outdir/$hemi.$mapname.sm00.$FSF_OUTPUT_FORMAT
        set cmd = (isxconcat0-sess -o $outvol $sessargs \
          -surf-anat $mapname -hemi $hemi)
        if($#flac == 0) set cmd = ($cmd -a $analysis)
        if($#flac)      set cmd = ($cmd -flac $flac -fsd $fsd)
        if($#tmpdir)  set cmd = ($cmd -tmp $tmpdir);
        if($paireddiff) set cmd = ($cmd -paired-diff)
        if($pairedavg)  set cmd = ($cmd -paired-avg)
        echo "cd `pwd`" |& tee -a $LF
        echo $cmd       |& tee -a $LF
        $cmd            |& tee -a $LF
        if($status) exit 1;
      end # Loop over maps
    else
      set outvol = $outdir/anat.norm.sm00.$FSF_OUTPUT_FORMAT
      set cmd = (isxconcat0-sess -o $outvol $sessargs \
         -vol-anat norm)
      if($#flac == 0) set cmd = ($cmd -a $analysis)
      if($#flac)      set cmd = ($cmd -flac $flac -fsd $fsd)
      if($#tmpdir)  set cmd = ($cmd -tmp $tmpdir);
      if($paireddiff) set cmd = ($cmd -paired-diff)
      if($pairedavg)  set cmd = ($cmd -paired-avg)
      echo "cd `pwd`" |& tee -a $LF
      echo $cmd       |& tee -a $LF
      $cmd            |& tee -a $LF
      if($status) exit 1;
    endif
  end # Loop over hemis
   
endif

echo "" | tee -a $LF
echo "" | tee -a $LF
echo "Started at $StartTime" | tee -a $LF
echo "Ended   at `date`"     | tee -a $LF
echo "isxconcat-sess completed" | tee -a $LF
echo " "
echo " "

exit 0;
###############################################

############--------------##################
parse_args:
set cmdline = ($argv);
while( $#argv != 0 )

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "-o":
      if($#argv == 0) goto arg1err;
      set outtop = $argv[1]; shift;
      breaksw

    case "-a":
    case "-analysis":
      if ( $#argv == 0) goto arg1err;
      set analysis = $argv[1]; shift;
      breaksw

    case "-f":
    case "-flac":
      if ( $#argv == 0) goto arg1err;
      set flac = $argv[1]; shift;
      breaksw

    case "-fsd":
      if ( $#argv == 0) goto arg1err;
      set fsd = $argv[1]; shift;
      breaksw

    case "-c":
    case "-contrast":
      if ( $#argv == 0) goto arg1err;
      set contrastlist = ($contrastlist $argv[1]); shift;
      breaksw

    case "-all-contrasts":
    case "-call":
      set DoAllContrasts = 1;
      breaksw

    case "-pct"
    case "-percent"
      set DoPercent = 1;
      breaksw

    case "-m":
    case "-map":
      if ( $#argv == 0) goto arg1err;
      set conmaplist = ($conmaplist $argv[1]); shift;
      breaksw

    case "-interp":
      if ( $#argv == 0) goto arg1err;
      set interp = $argv[1]; shift;
      breaksw

    case "-reg":
      if ( $#argv == 0) goto arg1err;
      set regfile = $argv[1]; shift;
      breaksw

    case "-fwhm":
      if ( $#argv == 0) goto arg1err;
      set fwhm = $argv[1]; shift;
      breaksw

    case "-talres":
      if ( $#argv == 0) goto arg1err;
      set talres = $argv[1]; shift;
      set ToVol = 1;
      breaksw

    case "-talxfm":
      if ( $#argv == 0) goto arg1err;
      set talxfm = $argv[1]; shift;
      set ToVol = 1;
      breaksw

    case "-morph":
      set DoMorph = 1;
      set ToVol = 1;
      breaksw

    case "-hemis":
      set hemilist = (lh rh); 
      set ToSurf = 1;
      breaksw

    case "-hemi":
    case "-h":
      if ( $#argv == 0) goto arg1err;
      set hemilist = ($hemilist $argv[1]); shift;
      set ToSurf = 1;
      breaksw

    case "-projfrac":
      if ( $#argv == 0) goto arg1err;
      set projfrac = $argv[1]; shift;
      set ToSurf = 1;
      breaksw

    case "-projfrac-avg":
      set DoProjFracAvg = 1;
      set ToSurf = 1;
      breaksw

    case "-trgsubject":
    case "-subject":
      if( $#argv == 0) goto arg1err;
      set trgsubject = $argv[1]; shift;
      set ToSurf = 1;
      breaksw

    case "-fsgd":
      if($#argv == 0) goto arg1err;
      set fsgd = $argv[1]; shift;
      if(! -e $fsgd) then
        echo "ERROR: cannot find $fsgd"
        exit 0;
      endif
      breaksw

    case "-tmpdir":
      if ( $#argv == 0) goto arg1err;
      set tmpdir = $argv[1]; shift;
      set cleanup = 0;
      breaksw

    case "-paired-diff"
      set paireddiff = 1;
      breaksw

    case "-paired-avg"
      set pairedavg = 1;
      breaksw

    case "-mean"
      set mean = 1;
      breaksw

    case "-no-anamaps"
      set GetHRF = 0;
      set anamaplist = ();
      breaksw

    case "-no-hdr"
    case "-no-hrf"
    case "-no-fir"
      set GetHRF = 0;
      breaksw

    case "-fir"
    case "-hdr"
    case "-hrf"
      set GetHRF = 1;
      breaksw

    case "-anat"
    case "-anats"
      set DoAnats = 1;
      breaksw

    case "-no-anats"
      set DoAnats = 0;
      breaksw

    case "-osgm"
      set DoOSGM = 1;
      set DoWRFX = 1;
      set DoRFX = 1;
      set DoFFX = 1;
      breaksw

    case "-wrfx"
      set DoWRFX = 1;
      set DoOSGM = 1;
      breaksw

    case "-rfx"
      set DoRFX = 1;
      set DoOSGM = 1;
      breaksw

    case "-ffx"
      set DoFFX = 1;
      set DoOSGM = 1;
      breaksw

    case "-split-mvc"
      set SplitMVC = 1;
      breaksw

    case "-no-split-mvc"
      set SplitMVC = 0;
      breaksw

    case "-debug":
      set verbose = 1;
      set echo = 1;
      breaksw

    case "-cwd":
      breaksw

    case "-s":
    case "-sf":
    case "-df":
    case "-d":
      set sessargs = ($sessargs $flag $argv[1]); shift;
      breaksw

    default:
      echo ERROR: Flag $flag unrecognized. 
      echo $cmdline
      exit 1
      breaksw
  endsw

end

goto parse_args_return;
############--------------##################

############--------------##################
check_params:

  if($#analysis == 0 && $#flac == 0) then
    echo "ERROR: no analysis specified"
    exit 1;
  endif

  if($#analysis != 0 && $#flac != 0) then
    echo "ERROR: cannot specify flac and analysis"
    exit 1;
  endif

  if($#analysis) then
    if(! -e $analysis) then
      echo "ERROR: cannot find $analysis"
      exit 1;
    endif
    set IsNative = `getana -a $analysis -t IsNative`
    set fsd = `getana -a $analysis -t fsd`;
    set IsSurface = `getana -a $analysis -t IsSurface`
    if($IsSurface) then
      set ToSurf = 1;
      set ToVol = 0;
      set trgsubject = `getana -a $analysis -t subject`
      if($trgsubject == self) then 
        echo "ERROR: analysis space is self surface, not supported"
        exit 1;
      endif
      set hemilist  = `getana -a $analysis -t hemi`
      set DoResample = 0;
    endif
    set UseTalairach = `getana -a $analysis -t IsMNI305`;
    if($UseTalairach) then
      set ToSurf = 0;
      set ToVol = 1;
      set hemilist = (mni305)
      set DoResample = 0;
    endif
    # Get the h-file, usually for FIRs
    if($#GetHRF == 0) set GetHRF = `getana -a $analysis -t IsFIR`;
    if($GetHRF) set anamaplist = ($anamaplist h);
  endif
  #echo "DoResample $DoResample"

  if($#flac) then
    set analysis = `basename $flac .flac`
    if($#fsd == 0) then
      echo "ERROR: need fsd with flac" 
      exit 1;
    endif
  endif

  if($DoAllContrasts) then
    # Wont include merged (ie, conjunction) contrasts
    set contrastlist = `getana -a $analysis -t contrasts`
  endif
  if($#contrastlist == 0) then
    echo "ERROR: no contrast specified"
    exit 1;
  endif
  if($#conmaplist == 0) then
    if(! $DoPercent ) then
      set conmaplist = (ces cesvar)
    else
      set conmaplist = (cespct cesvarpct)
    endif
  endif

  if($#outtop == 0) then
    echo "ERROR: no output directory specified"
    exit 1;
  endif
  mkdir -p $outtop
  if($status) then
    echo "ERROR: making $outtop";
    exit 1;
  endif

  if($paireddiff && $pairedavg) then
    echo "ERROR: cannot --paired-diff and --paired-avg";
    exit 1;
  endif

  if($?FSF_OUTPUT_FORMAT == 0) setenv FSF_OUTPUT_FORMAT mgh

  if($DoResample) then
    if($ToSurf == 0 && $ToVol == 0) set ToVol = 1;
    if($ToSurf && $ToVol) then
      echo "ERROR: cannot convert to surface and volume"
      exit 1;
    endif
    if("$interp" != "trilin" && "$interp" != "nearest") then
      echo "ERROR: interp is $interp, must be trilinear or nearest"
      exit 1;
    endif
    if($talres != 1 && $talres != 2) then
      echo "ERROR: talres = $talres, must be 1 or 2"
      exit 1;
    endif
    if($#hemilist) then
      foreach hemi ($hemilist)
        if($hemi != lh && $hemi != rh) then
          echo "ERROR: hemi = $hemi, must be rh or lh"
          exit 1;
        endif
      end
    else
      set hemilist = mni305;
    endif
    if($#projfrac && $DoProjFracAvg) then
      echo "ERROR: cannot --projfrac and --projfrac-avg"
      exit 1
    endif
    # default proj frac
    if($#projfrac == 0 && $DoProjFracAvg == 0) then
      set projfrac = 0.5;
    endif
  endif
  if($#regfile == 0) set regfile = register.dof$RegDOF.dat 

  if($talres == 2) then
    set subcortmask = $SUBJECTS_DIR/fsaverage/mri.2mm/subcort.mask.mgz
  else
    set subcortmask = $SUBJECTS_DIR/fsaverage/mri/subcort.mask.1mm.mgz
  endif
  if(! -e $subcortmask) then
    echo "ERROR: cannot find $subcortmask"
    exit 1;
  endif


goto check_params_return;
############--------------##################

############--------------##################
arg1err:
  echo "ERROR: flag $flag requires one argument"
  exit 1
############--------------##################

############--------------##################
usage_exit:
  echo " "
  echo "USAGE: isxconcat-sess"
  echo ""
  echo "   -o outdir       : top group output directory"
  echo ""
  echo "   -a analysisname : session-level functional analysis name"
  echo "   -c contrast     : contrast "
  echo "   -all-contrasts  : contrast "
  echo ""
  echo "   -paired-diff    : compute 1-2, 3-4, etc"
  echo "   -no-hrf         : create hrf maps (saved as space.h.nii)"
  echo "   -no-anamaps     : do not compute analysis level maps"
  echo "   -no-split-mvc   : do not split multi-variate contrasts"
  echo ""
  echo "   -percent        : use voxel-wise percent signal change"
  echo "   -m mapname      : contrast-level maps stem (def is ces and cesvar)"
  echo ""
  echo "   -paired-avg     : compute avg 1+2, 3+4, etc"
  echo "   -mean           : mean over all frames"
  echo ""
  echo " Session arguments"
  echo "   -sf sessidfile  ..."
  echo "   -df srchdirfile ..."
  echo "   -s  sessid      ..."
  echo "   -d  srchdir     ..."
  echo ""
  echo "   -fsgd fsgdfile : for comparing subject order"
  echo ""
  echo "   -help"
  echo "   -debug"
  echo "   -tmpdir dir : use dir for temporary storage (implies -nocleanup)"
  echo "   -nocleanup  : do not delete tmpdir"
  echo "   -version    : print version and exit"
  echo " "
  echo "Options below are for analyses done in native space"
  echo "   -reg regfile    : instead of register.dat"
  echo "   -interp method  : trilinear (def) or nearest"
  echo "   -fwhm fwhm      : smooth by fwhm"
  echo ""
  echo "  Volume-based options:"
  echo "   -talres res     : 1 or 2 mm (def is 2)"
  echo "   -talxfm xfmfile : use xfmfile instead of talairach.xfm"
  echo ""
  echo "  Surface-based options:"
  echo "   -hemis           : surface-based, both lh and rh";
  echo "   -hemi hemisphere : for surface-based analysis (lh or rh)";
  echo "   -projfrac frac   : sample frac into surface bet white and gray (0.5 def)"
  echo "   -projfrac-avg    : average across cortical depth"
  echo "   -trgsubject subject : sample to subject's surface (def fsaverage)"
  echo ""

  if(! $PrintHelp ) exit 1;

  echo " "
  echo "$VERSION "
  echo " "

  cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
exit 1;

#---- Everything below here is printed out as part of help -----#
BEGINHELP

Creates the output necessary to run group fMRI analysis by
concatenating the results from a first-level analysis of all the
sessions together into one multi-frame file suitable for use with
mri_glmfit. It also creates several other files that can be useful for
analysis and quality control. The output is stored in a directory
structure that closely follows the FS-FAST directory structure.

-o outdir

Output directory. The data will be placed in outdir/analysis. All
analyses that are placed in the directory should have the same
set of sessions in the same order.

-a analysis

First-level analysis made with mkanalysis-sess (and analyzed with
selxavg3-sess).

-c contrast

First-level contrast made with mkcontrast-sess. You can specify
multiple contrasts with multiple -c arguments. Or use -all-contrasts.

-paired-diff

Compute the difference between adjacent sessions as a preprocessing
step to performing a paired-t analysis. This computes the differences
between the first session and the second, the third and the forth, etc.
The sign first-second, eg: 1-2 3-4 5-6. This reduces the number of 
output frames by a factor of 2. See also -paried-avg.

-no-hrf

For FIR designs, do not create an average hemodynamic response file.
This file can be quite large and take a long time to produce.

-no-anamaps

Only computes maps based on the contrast. By default, there are 
several maps that are computed that are based on just the 
analysis (eg, meanfunc is the mean functional). 

-no-split-mvc

Do not split multi-variate contrasts (MVCs). By default, each variate
is split into a separate file (eg, ces.000.nii, ces.001.nii, etc).
This causes them to simply be concatenated together.

-percent

Use cespct instead of ces. This is the percent contrast computed 
on a voxel-wise basis.

-m mapname

Use mapname as the map in the contrast. Default is ces. This is
mostly for experimental purposes.

-paired-avg

Compute the average adjacent sessions as a preprocessing step. This
computes the average of the first session and the second, the third
and the forth, etc.  eg: (1+2)/2 (3+4)/2 (5+6)/2. This reduces the
number of output frames by a factor of 2. See also -paried-diff.

OUTPUT 

All the data will be under a folder called outdir/analysis. This 
directory will have a subfolder for each contrast. These
are some of the output files and what they mean:

sessid.txt - list of the sessions specified

subjectlist.txt - list of subjects specified. This can be used with
 mris_preproc, asegstats2table, aparcstats2table, or
 make_average_subject.

analysis.info - copy of the file created by mkanalysis-sess

ffxdof.dat - total number of degrees of freedom to use for a fixed
  effects analysis

masks.nii - maps of the binary mask for each subject

mask.nii - intersection of the binary masks for each subject

meanfunc.nii - mean functional intensity maps for each subject

rho1.nii - temporal AR1 for each subject

raw.fsnr.nii - maps of the raw functional SNR for each subject

fsnr.nii - maps of the functional SNR for each subject computed after processing

CONTRAST OUTPUT

There will be a subfolder for each contrast. There will be several
files in this folder:

contrast.lowerlevel.mat -- this is a copy of the contrast mat file
  created by mkcontrast-sess.

ces.nii - "contrast effect size". This is maps of contrast values,
  one map for eacah subject. This is the file that will be used
  for group analysis. 

cesvar.nii - This is maps of the variance of contrast values,
  one map for eacah subject. This is the file that will be used
  for weighted random and fixed effects group analysis.