#! /bin/tcsh -f

#
# pctsurfcon
#
# Original Author: Doug Greve
#
# Copyright © 2021 The General Hospital Corporation (Boston, MA) "MGH"
#
# Terms and conditions for use, reproduction, distribution and contribution
# are found in the 'FreeSurfer Software License Agreement' contained
# in the file 'LICENSE' found in the FreeSurfer distribution, and here:
#
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
#
# Reporting: freesurfer@nmr.mgh.harvard.edu
#
#

#
set VERSION = 'pctsurfcon dev';
set inputargs = ($argv);

set subject = ();
set fsvol = rawavg;
set mov = ();
set reg = ();
set UseMask = 1;

set DoLH = 1;
set DoRH = 1;

set Reshape = (--reshape);
set Reshape = (--noreshape);

set GMProjFrac = ();
set GMProjAbs = ();
set WMProjAbs = 1;

set Interp = trilinear
set OutBase = w-g.pct;
set ConSign = 100;
set surf = ()
set ext = mgh;

set debug = 0;
set tmpdir = ();
set cleanup = 1;
set PrintHelp = 0;
set nolog = 0;

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

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

set outdir = $SUBJECTS_DIR/$subject/surf
if($#tmpdir == 0) set tmpdir = $outdir/tmp.pctsurfcon.$$
mkdir -p $tmpdir

if(! $nolog) then
  set LF = $SUBJECTS_DIR/$subject/scripts/pctsurfcon.log
  if(-e $LF) mv $LF $LF.old
else
  set LF = /dev/null
endif

echo "Log file is $LF"

echo "Logfile for pctsurfcon" >> $LF
date |& tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" |& tee -a $LF
echo "cd `pwd`" |& tee -a $LF
echo $0 |& tee -a $LF
echo $VERSION |& tee -a $LF
uname -a |& tee -a $LF
echo "FREESURFER_HOME $FREESURFER_HOME" |& tee -a $LF

set hemilist = ();
if($DoLH) set hemilist = ($hemilist lh);
if($DoRH) set hemilist = ($hemilist rh);

foreach hemi ($hemilist)

  # Sample surf in white matter (inside the surface)
  set wm = $tmpdir/$hemi.wm.$ext
  set cmd = (mri_vol2surf --mov $mov \
    --hemi $hemi $Reshape --interp $Interp \
    --projdist -$WMProjAbs --o $wm)
  if($#reg == 0) set cmd = ($cmd --regheader $subject);
  if($#reg != 0) set cmd = ($cmd --reg $reg);
  if($UseMask) set cmd = ($cmd --cortex);
  if($#surf) set cmd = ($cmd --surf $surf)
  echo $cmd
  $cmd
  if($status) exit 1;

  # Sample surf in gray matter (outside the surface)
  set gm = $tmpdir/$hemi.gm.$ext
  set cmd = (mri_vol2surf --mov $mov \
    --hemi $hemi $Reshape --interp $Interp --o $gm)
  if($#GMProjFrac) set cmd = ($cmd --projfrac $GMProjFrac)
  if($#GMProjAbs) set cmd = ($cmd --projdist $GMProjAbs)
  if($#reg == 0) set cmd = ($cmd --regheader $subject);
  if($#reg != 0) set cmd = ($cmd --reg $reg);
  if($UseMask) set cmd = ($cmd --cortex);
  if($#surf) set cmd = ($cmd --surf $surf)
  echo $cmd
  $cmd
  if($status) exit 1;

  # Compute percent difference pct = 100*(w-g)/((w+g)/2)
  set out = $outdir/$hemi.$OutBase.$ext
  set cmd = (mri_concat $wm $gm --paired-diff-norm \
    --mul $ConSign --o $out)
  echo $cmd 
  $cmd
  if($status) exit 1;

  # Compute stats
  set sum = $SUBJECTS_DIR/$subject/stats/$hemi.$OutBase.stats
  set cmd = (mri_segstats --in $out --annot $subject $hemi aparc \
    --sum $sum --snr);
  echo $cmd 
  $cmd
  if($status) exit 1;

end

# Cleanup
if($cleanup) then
  echo "Cleaning up" |& tee -a $LF
  rm -r $tmpdir
endif

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

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

  set flag = $argv[1]; shift;

  switch($flag)

    case "--s":
      if ( $#argv < 1) goto arg1err;
      set subject = $argv[1]; shift;
      set subject = `basename $subject`; # remove trailing "/"
      breaksw

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

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

    case "--reg":
      if ( $#argv < 1) goto arg1err;
      set reg = $argv[1]; shift;
      if(! -e $reg) then
        echo "ERROR: cannot find $reg"
        exit 1;
      endif
      #set subject = `head -n 1 $reg`
      set subject = `reg2subject --r $reg`;
      breaksw

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

    case "--no-mask":
      set UseMask = 0;
      breaksw

    case "--lh-only":
    case "--no-rh":
      set DoRH = 0;
      breaksw

    case "--rh-only":
    case "--no-lh":
      set DoLH = 0;
      breaksw

    case "--gm-proj-frac":
      if($#argv < 1) goto arg1err;
      set GMProjFrac = $argv[1]; shift;
      breaksw

    case "--gm-proj-abs":
      if($#argv < 1) goto arg1err;
      set GMProjAbs = $argv[1]; shift;
      breaksw

    case "--wm-proj-abs":
      if($#argv < 1) goto arg1err;
      set WMProjAbs = $argv[1]; shift;
      breaksw

    case "--mgh":
      set ext = mgh;
      breaksw
    case "--mgz":
      set ext = mgz;
      breaksw
    case "--nii.gz":
      set ext = nii.gz;
      breaksw

    case "--nearest":
      set Interp = nearest
      breaksw

    case "--trilin":
    case "--trilinear":
      set Interp = trilinear
      breaksw

    case "--neg":
      set ConSign = -100;
      breaksw

    case "--pos":
      set ConSign = 100;
      breaksw

    case "--pial":
      set surf = pial
      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 "--debug":
      set verbose = 1;
      set echo = 1;
      breaksw

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

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

end

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

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

  if($#subject == 0) then
    echo "ERROR: must spec a subject id"
    exit 1;
  endif
  if(! -e $SUBJECTS_DIR/$subject) then
    echo "ERROR: cannot find $subject in $SUBJECTS_DIR"
    exit 1;
  endif

  if($#mov == 0) set mov = $SUBJECTS_DIR/$subject/mri/$fsvol.mgz
  if(! -e $mov) then
    echo "ERROR: cannot find $mov"
    exit 1;
  endif

  if($#GMProjFrac && $#GMProjAbs) then
    echo "ERROR: cannot --gm-proj-abs and --gm-proj-frac"
    exit 1;
  endif

  if($#GMProjFrac == 0 && $#GMProjAbs == 0) set GMProjFrac = 0.3;

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 "pctsurfcon - compute surface-wise gray/white contrast"
  echo ""
  echo "Required Arguments:";
  echo ""
  echo "  --s subject  : FreeSurfer subject name"
  echo ""
  echo "Other Arguments:";
  echo ""
  echo "  --fsvol fsvol : use fsvol instead of $fsvol"
  echo "  --b outbase : use outbase instead of $OutBase (?h.$OutBase.mgh)"
  echo "  --lh-only, --no-rh : compute lh only"
  echo "  --rh-only, --no-lh : compute rh only"
  echo "  --gm-proj-frac frac : GM projection fraction (default, 0.3)"
  echo "  --gm-proj-abs  dist : GM projection distance (default is to use frac)"
  echo "  --wm-proj-abs  dist : WM projection distance (default is $WMProjAbs mm)"
  echo "  --neg : compute G-W instead of W-G"
  echo "  --no-mask : do not mask out non-cortical regions"
  echo "  --pial : use pial surface as base (instead of white) to compute gray/CSF contrast"
  echo ""
  echo "Even more arguments"
  echo ""
  echo " --mgz : extension, default ext is $ext"
  echo " --nii.gz : extension, default ext is $ext"
  echo ""
  echo "  --tmp tmpdir  : temporary dir (implies --nocleanup)"
  echo "  --nocleanup  : do not delete temporary files"
  echo "  --version : print version and exit"
  echo "  --help    : print help and exit"
  echo ""

  if($PrintHelp) \
  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

Computes the vertex-by-vertex percent contrast between white and gray
matter. The computation is:

         100*(W-G)
   pct = ---------
         0.5*(W+G)

The white matter is sampled 1mm below the white surface. This distance
is changeable with --wm-proj-abs. The gray matter is sampled 30% the
thickness into the cortex. This fraction is changeable with
--gm-proj-frac. The volume that is sampled is rawavg.mgz, but this
can be changed with --fsvol.

The output is stored in the surf dir of the given subject as
?h.w-g.pct.mgh. This can be changed to ?h.outbase.mgh with 
--o outbase.

Note: the non-cortex region is zeroed-out by default. This can be
overridden with the --no-mask flag.

Example:

pctsurfcon --s subject