#!/bin/tcsh -f
# groupstatsdiff

set VERSION = 'groupstatsdiff 8.2.0';

set g1dir = ();
set g2dir = ();
set DoMaps = 1;
set tmpdir = ();
set cleanup = 1;
set LF = ();
set outdir = ()
set surflist = (white)
set DoOSGM = 0;
set DoCommon = 1;
set AllowDiffSubjs = 0; # allow subject list to be diff
set nomeaslist = ();
set NoBA = 0;
set DoAparcStats = 1
set DoAsegStats = 1
set DoWMParcStats = 1
set DoPrune = 1
set fwhmlist = ()
set sd1 = ()
set sd2 = ()
set dicectab = $FREESURFER_HOME/ASegStatsLUT.txt
set DoDice = 1;

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

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

set StartTime = `date`;
set tSecStart = `date '+%s'`;

mkdir -p $outdir/rois $outdir/maps
pushd $outdir > /dev/null
set outdir = `pwd`;
popd > /dev/null

# Set up log file
if($#LF == 0) set LF = $outdir/groupstatsdiff.log
if($LF != /dev/null) rm -f $LF
echo "Log file for groupstatsdiff" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
echo "" | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
uname -a  | tee -a $LF

set mapmeaslist  = `cat $g1dir/mapmeas.list.txt`
set hemilist     = `cat $g1dir/hemi.list.txt`
set aparclist    = `cat $g1dir/aparc.list.txt`
if($#fwhmlist == 0) set fwhmlist     = `cat $g1dir/fwhm.list.txt`

if($DoOSGM) then
   set gcmd = (--osgm)
else
  set gcmd = (--fsgd $fsgd)
endif

#========================================================

# ASEG Stats ----------------------------------------------------
if($DoAsegStats) then
foreach hemi ($hemilist other)
  foreach meas (volume mean)
    set b = aseg.$hemi.$meas
    if($meas == mean) set b = aseg.$hemi.intensity
    set table1 =  $g1dir/rois/$b.dat
    set table2 =  $g2dir/rois/$b.dat
    set difftable = $outdir/rois/diff.$b.dat
    set cmd = (stattablediff --t1 $table1 --t2 $table2 --o $difftable --percent)
    if($DoCommon)    set cmd = ($cmd --common)
    if($AllowDiffSubjs) set cmd = ($cmd --diff-subjs) # allow subject list not the same in both
    echo "#C# $cmd" | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
    set glmdir = $outdir/rois/glm.$b
    set cmd = (mri_glmfit --table $difftable --glmdir $glmdir $gcmd)
    echo "#C# $cmd"  | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
  end # meas
end # hemi

if($DoDice) then
  set slist = (`grep Input $fsgd| awk '{print $2}'`)
  set dicedir = $outdir/rois/aseg.dice
  mkdir -p $dicedir
  set dicefile = $dicedir/dice.dat 
  cp $dicectab $dicedir/aseg.ctab
  rm -f $dicefile $dicedir/slist.txt
  foreach s ($slist)
    echo $s >> $dicedir/slist.txt
    set aseg1 = $sd1/$s/mri/aseg.mgz
    set aseg2 = $sd2/$s/mri/aseg.mgz
    set dicedat   = $dicedir/$s.dice.dat
    set dicetable = $dicedir/$s.dice.table.dat
    set cmd = (mri_compute_seg_overlap -dice $aseg1 $aseg2 $dicectab 1 0 $dicedat $dicetable)
    echo "$cmd" | tee -a $LF
    $cmd |& tee -a $LF 
    if($status) goto exit 1;
    cat $dicedat | tee -a $dicefile
  end
endif # DoDice
endif # DoAseg

# APARC Stats ----------------------------------------------------
if($DoAparcStats) then
foreach aparc ($aparclist)
  foreach surf ($surflist)
    if($aparc == aparc.a2009s && $surf == pial) continue;
    if($NoBA  && $aparc == BA) continue;
    foreach hemi ($hemilist)
      foreach meas ($mapmeaslist)
        if($meas == curv || $meas == sulc) continue # curv breaks aparcstats2table
        if($meas == w-g.pct && ($aparc != aparc || $surf != white)) continue 
        set Skip = 0;
        foreach nomeas ($nomeaslist)
          if($meas == $nomeas) then
            set Skip = 1;
            break; 
          endif
        end
	if($Skip) continue;
        set parc = $aparc
        if($surf == pial) set parc = pial.$aparc; # only aparc
        set table1 = $g1dir/rois/$parc.$hemi.$surf.$meas.dat
        set table2 = $g2dir/rois/$parc.$hemi.$surf.$meas.dat
        set difftable = $outdir/rois/diff.$parc.$hemi.$surf.$meas.dat
        # --rm-exvivo removes the string 'exvivo' from the ROI name
        # to make 6.0 compatible with 5.3
        set cmd = (stattablediff --t1 $table1 --t2 $table2 --o $difftable --percent)
        set cmd = ($cmd --rm-exvivo)
        if($DoCommon) set cmd = ($cmd --common)
        if($AllowDiffSubjs) set cmd = ($cmd --diff-subjs) # allow subject list not the same in both
        echo "#C# $cmd" | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
        set glmdir = $outdir/rois/glm.$parc.$hemi.$surf.$meas
        set cmd = (mri_glmfit --table $difftable --glmdir $glmdir $gcmd)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
      end # meas
    end # hemi
  end # surf
end # aparc
endif


# WMParc Stats --------------------------------------------
if($DoWMParcStats) then
set table1 =  $g1dir/rois/wmparc.vol.dat
set table2 =  $g2dir/rois/wmparc.vol.dat
set difftable = $outdir/rois/diff.wmparc.vol
set cmd = (stattablediff --t1 $table1 --t2 $table2 --o $difftable --percent)
if($DoCommon) set cmd = ($cmd --common)
if($AllowDiffSubjs) set cmd = ($cmd --diff-subjs) # allow subject list not the same in both
echo "#C# $cmd" | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
set glmdir = $outdir/rois/glm.wmparc.vol
set cmd = (mri_glmfit --table $difftable --glmdir $glmdir $gcmd)
echo "#C# $cmd"  | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
endif

# Maps and GLM ---------------------------------------
if($DoMaps) then
  foreach meas ($mapmeaslist)
  set Skip = 0;
    foreach nomeas ($nomeaslist)
      if($meas == $nomeas) then
        set Skip = 1;
        break;
      endif
    end
    if($Skip) continue;
    foreach hemi ($hemilist)
      foreach fwhm ($fwhmlist)
        set fwhmstr = `printf %02d $fwhm`
        set sm1file =  $g1dir/maps/$meas.$hemi.sm$fwhmstr.mgh
        set sm2file =  $g2dir/maps/$meas.$hemi.sm$fwhmstr.mgh
        # Compute the difference, set to 0 any voxel where either input is 0
        set smdifffile = $outdir/maps/diff.$meas.$hemi.sm$fwhmstr.mgh	
        set cmd = (fscalc $sm1file sub0 $sm2file -o $smdifffile)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
        # glm ------------------------
        set glmdir = $outdir/maps/glm.$meas.$hemi.sm$fwhmstr
        set cmd = (mri_glmfit --surface fsaverage $hemi --y $smdifffile --glmdir $glmdir $gcmd)
        if($DoPrune == 0) set cmd = ($cmd --no-prune)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) then
	  if($meas == area) then
            # This can happen when comparing t1-only with the t2 stream
            echo "Though there was an error, it is probably due to pruning of area diffs" | tee -a $LF
            echo "so I'm just going to move on to the next one" | tee -a $LF
	    continue;
          endif
          exit 1;
        endif
	if($fwhm != 0) then
          set cmd = (mri_glmfit-sim --glmdir $glmdir --cwp .05 --cache 1.3 abs --2spaces)
          echo "#C# $cmd"  | tee -a $LF
          $cmd | tee -a $LF
          if($status) exit 1;
        endif
        # Compute the percent difference, set to 0 any voxel where either input is 0
        set smdifffile = $outdir/maps/pctdiff.$meas.$hemi.sm$fwhmstr.mgh	
        set cmd = (fscalc $sm1file pctdiff0 $sm2file -o $smdifffile)
        echo "#C# $cmd"  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
      end # fwhm
    end # hemi
  end #meas
endif

#========================================================
# Done

echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "groupstatsdiff-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "groupstatsdiff Done" |& tee -a $LF

exit 0

###############################################

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

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

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

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

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

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

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

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

    case "--no-dice":
      set DoDice = 0;
      breaksw

    case "--osgm":
      set DoOSGM = 1;
      breaksw

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

    case "--common":
      set DoCommon = 1;
      breaksw
    case "--no-common":
      set DoCommon = 0;
      breaksw

    case "--allow-subj-diff":
      set AllowDiffSubjs = 1;
      breaksw
    case "--no-allow-subj-diff":
      set AllowDiffSubjs = 0;
      breaksw

    case "--no-maps":
      set DoMaps = 0;
      breaksw

    case "--no-area":
    case "--noarea":
      set nomeaslist = ($nomeaslist area);
      breaksw

    case "--no-volume":
    case "--novolume":
      set nomeaslist = ($nomeaslist volume)
      breaksw

    case "--no-aparcstats":
      set DoAparcStats = 0
      breaksw

    case "--no-asegstats":
      set DoAsegStats = 0
      breaksw

    case "--no-wmparcstats":
      set DoAsegStats = 0
      breaksw

    case "--no-stats":
      set DoAparcStats = 0
      set DoAsegStats = 0
      set DoWMParcStats = 0
      breaksw

    case "--no-BA":
    case "--no-ba":
    case "--noba":
      set NoBA = 1;
      breaksw

    case "--no-prune":
      set DoPrune = 0
      breaksw

    case "--lh":
      set hemilist = (lh)
      breaksw

    case "--rh":
      set hemilist = (rh)
      breaksw

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

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      breaksw

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

    case "--nocleanup":
      set cleanup = 0;
      breaksw

    case "--cleanup":
      set cleanup = 1;
      breaksw

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

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

end

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

############--------------##################
check_params:
if($#g1dir == 0) then
  echo "ERROR: must spec --g1"
  exit 1;
endif
if($#g2dir == 0) then
  echo "ERROR: must spec --g2"
  exit 1;
endif
if($#outdir == 0) then
  echo "ERROR: Must specify output folder"
  exit 1;
endif
set fsgd1 = $g1dir/group.fsgd
set fsgd2 = $g2dir/group.fsgd
set n = `diff $fsgd1 $fsgd2 | wc -l`
if(0 && $n != 0) then
  echo "ERROR: fsgd files $fsgd1 and $fsgd2 are not the same"
  exit 1; 
endif
set fsgd = $g1dir/group.fsgd

if($#sd1 == 0) set sd1 = `dirname $g1dir`
if($#sd2 == 0) set sd2 = `dirname $g2dir`

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

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

############--------------##################
usage_exit:
  echo ""
  echo "groupstatsdiff "
  echo "  --g1 group1dir : output folder from groupstats"
  echo "  --g2 group2dir : output folder from groupstats"
  echo "  --o output : output folder for difference"
  echo "  --no-maps : only analyze ROI data"
  echo "  --osgm : use OSGM instead of native FSGD"
  echo "  --no-common : do not select common segs when running stattablediff"
  echo "  --allow-subj-diff : allow list of subjects to be different between the two"
  echo "  --no-area : do not compute area diffs"
  echo "  --no-volume : do not compute volume diffs"
  echo "  --no-ba : do not compute diffs for BA labels"
  echo "  --no-aparcstats : do not do aparcstats"
  echo "  --no-asegstats : do not do asegstats"
  echo "  --no-wparcstats : do not do wmparcstats"
  echo "  --no-stats : do not do any ROI stats"
  echo "  --no-prune : do not prune when running mri_glmfit (may be nec if g1=g2)"
  echo "  --fwhm fwhm <--fwhm fwhm> : override the fwhm from group ana (must be present)"
  echo "  --sd1 sd1, --sd2 sd2 : subjects dirs for computing dice (default is parent dir of groupdir)"
  echo "  --no-dice : do not compute dice"
  echo "  --dice-ctab ctab : ctab to use for dice (default is $dicectab)"
  echo "  "

  if(! $PrintHelp) exit 1;
  echo $VERSION
  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

This is a script run on two outputs of groupstats with the intention
of evaluating the differences between two different recon-all analyses
(eg, two versions of FS, same version but different platforms, same
version but different parameters). The groupstats script is run
separately on each of the analyses, storing the output in group1dir
and group2dir which are then passed to this script.

This script can be run in two different modes, either with --osgm
(one-sample-group-mean) or without in which case it will use the FSGD
file from the groupstats run. If groupstats was run with OSGM, then
these two are the same. Use OSGM when simply looking at the difference
between the two analyses. Use the native FSGD when look at whether
there is an interaction between effect of interest (eg, age) and the 
two analyses.

One can then evaluate the differences with something like

tksurfer fsaverage lh inflated -aparc \
   -ov diffdir/maps/glm.thickness.lh.sm10/osgm/sig.mgh \
   -t diffdir/maps/diff.thickness.lh.sm10.mgh

The map will indicate OSGM differences and the "time" course will show
the paired differences in each subject. If you want to dig deeper into
a particular subject, then try using fvcompare.

Other maps that can be examined are: thickness, area, volume.