#! /bin/csh -f

#
# tkmedit-sess
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: greve $
#    $Date: 2014/04/04 18:27:27 $
#    $Revision: 1.21 $
#
# 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: tksurfer-sess,v 1.21 2014/04/04 18:27:27 greve Exp $'
set inputargs = ($argv);
set analysis      = ();
set contrastlist   = ();
set mapanalysis   = ();
set map = ();
set annot = aparc;
set hemi = ();
set surf = inflated
set ShowHDR = 1;
set PrintHelp = 0;
set fthresh = 2;
set fmax    = 4;
set DoAllContrasts = 0;
set TCL = ();
set FWHM = ();

set Delink = 0;
set UseVol2Surf = 1;
set ProjFrac = 0.5;
set Interp = nearest

set RawSpace = native;
set RawSpaceType = volume
set trgsubject = ();
set SurfReg = ();
set DoFieldSign = 0;
set UseTalairach = 0;
set RegDOF = 6;
set regfile = ();
set UseFreeview = 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;
  exit 1;
endif
set n = `echo $argv | grep -e -version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif
set SessList = (`getsesspath $argv`);
if($status || $#SessList == 0) then
  echo $SessList
  exit 1;
endif
goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

echo $VERSION
uname -a      
date          
which tksurfer

# ---------------------------------------------
foreach sess ($SessList)
  echo ""
  set sessid = `basename $sess`
  set self = `cat $sess/subjectname`;

  # Determine the source subject
  if($RawSpace == self || $RawSpace == native) then
    set srcsubject = $self
  else if($RawSpaceType == surface) then
    # Non-self surface space
    set srcsubject = $RawSpace
  else
    # Talairach Space
    set srcsubject = fsaverage
  endif
    
  # Determine target subject
  if($#trgsubject == 0) then
    if($RawSpaceType == surface) then
      if($RawSpace == self) then
        set trgsubjectuse = $self
      else
        set trgsubjectuse = $RawSpace
      endif
    else
      # Talairach Space
      set trgsubjectuse = fsaverage
    endif
  else
    # Target subject specified
    if($trgsubject == self) then
      set trgsubjectuse = $self
    else
      set trgsubjectuse = $trgsubject
    endif
  endif

  set anapath = $sess/$fsd/$analysis
  if(! -e $anapath) then
    echo "ERROR: cannot find $anapath"
    exit 1;
  endif

  # concatenate contrast files
  set mapanapath = $sess/$fsd/$mapanalysis
  set maplist = ()
  echo "list $contrastlist"
  foreach contrast ($contrastlist)
    set conpath = $mapanapath/$contrast
    if(! -e $conpath) then
      echo "ERROR: cannot find $conpath"
      exit 1;
    endif
    if($DoFieldSign && $RawSpace == native) then
      set mapstem = $conpath/$hemi.fieldsign
    else
      set mapstem = $conpath/$map
    endif
    set mappath = `stem2fname $mapstem`;
    if($status) then
      echo "$mappath"
      exit 1;
    endif
    set maplist = ($maplist $mappath)
    echo con/map $contrast $mappath
  end
  echo maplist $maplist
  if($?FREESURFER_TMP_DIR) then
    set mapconcatfile = $FREESURFER_TMP_DIR/$sessid.$mapfsd.$mapanalysis.tksconcat.nii
  else
    set mapconcatfile = $sess/$mapfsd/$mapanalysis/tksconcat.nii
  endif
  set cmd = (mri_concat $maplist --o $mapconcatfile);
  echo $cmd
  $cmd
  if($status) exit 1;

  # handle FIR
  if($ShowHDR) then
    set hstem = $anapath/h
    set hpath = `stem2fname $hstem`;
    if($status) then
      echo "$hpath"
      exit 1;
    endif
    set hoffsetstem = $anapath/h-offset
    set hoffsetpath = `stem2fname $hoffsetstem`;
    if($status) then
      echo "$hoffsetpath"
      exit 1;
    endif
  endif

  if($RawSpace == native) then 
    set reg = $sess/$fsd/$regfile
    if(! -e $reg) then
      echo "ERROR: cannot find $reg"
      exit 1
    endif
  endif

  if($RawSpaceType == volume && ! $DoFieldSign) then
    # Sample volume to the surface
    set reg = ();
    if($RawSpace == native) then 
      set reg = $sess/$fsd/$regfile
      if(! -e $reg) then
        echo "ERROR: cannot find $reg"
        exit 1
      endif
    endif
    # Set up vol2surf basic cmd
    set cmd0 = (mri_vol2surf --hemi $hemi --projfrac $ProjFrac \
      --interp $Interp --cortex)
    if($RawSpace == native)    set cmd0 = ($cmd0 --reg $reg)
    if($RawSpace == mni305) set cmd0 = ($cmd0 --regheader $srcsubject)
    if($#FWHM) set cmd0 = ($cmd0 --surf-fwhm $FWHM)
    # run vol2surf on the maps
    set mov = $mapconcatfile
    set out = $sess/$mapfsd/$mapanalysis/$hemi.tksconcat.mgh
    set cmd = ($cmd0 --mov $mov --o $out)
    echo $cmd
    $cmd
    if($status) exit 1;
    rm -f $mov
    set mapconcatfile = $out
    if($ShowHDR) then
      # HRF
      set mov = $hpath
      set out = $sess/$mapfsd/$mapanalysis/$hemi.h.mgh
      set cmd = ($cmd0 --mov $mov --o $out)
      echo $cmd
      $cmd
      if($status) exit 1;
      set hpath = $out
      # Offset
      set mov = $hoffsetpath
      set out = $sess/$mapfsd/$mapanalysis/$hemi.h-offset.mgh
      set cmd = ($cmd0 --mov $mov --o $out)
      if($RawSpace == native)    set cmd0 = ($cmd0 --reg $reg)
      if($RawSpace == mni305) set cmd0 = ($cmd0 --regheader $srcsubject)
      echo $cmd
      $cmd
      if($status) exit 1;
      set hoffsetpath = $out
      cp $sess/$mapfsd/$mapanalysis/h.dat $sess/$mapfsd/$mapanalysis/$hemi.h.dat
    endif
  endif

  if($srcsubject != $trgsubjectuse) then
    # Basic cmd line
    set cmd0 = (mri_surf2surf --srcsubject $srcsubject \
     --trgsubject $trgsubjectuse --hemi $hemi)
    if($#SurfReg) set cmd0 = ($cmd0 --surfreg $SurfReg)
    # Convert the maps
    set cmd = ($cmd0 --sval $mapconcatfile --tval $mapconcatfile)
    echo $cmd
    $cmd
    if($status) exit 1;
    if($ShowHDR) then
      # Convert the hrfs
      set cmd = ($cmd0 --sval $hpath --tval $hpath)
      echo $cmd
      $cmd
      if($status) exit 1;
      set cmd = ($cmd0 --sval $hoffsetpath --tval $hoffsetpath)
      echo $cmd
      $cmd
      if($status) exit 1;
    endif
    set subject = $self
  endif

  if($UseFreeview == 0) then
    set cmd = (tksurfer $trgsubjectuse $hemi $surf -overlay $mapconcatfile)
    set cmd = ($cmd -label-outline -colscalebarflag 1)
    if($ShowHDR) then
      set cmd = ($cmd -timecourse $hpath -timecourse-offset $hoffsetpath)
    endif
    if($#annot) set cmd = ($cmd -annot $annot)
    set maskopt = ()
    if($IsRet && ($map == angle || $map == fieldsign)) then
      set mapanadir = $sess/$mapfsd/$mapanalysis
      if($IsSurface) set maskstem = $mapanadir/fieldsign/fsig
      if(! $IsSurface) then
        set maskstem = $mapanadir/fieldsign/$hemi.fsig
        if($status) then
          echo "ERROR: cannot find $mapanadir/fieldsign/$hemi.fsig."
          echo "       Run fieldsign-sess."
          exit 1;
        endif
      endif
      set mask = `stem2fname $maskstem`
      set maskopt = (-mask $mask -mask-thresh $fthresh -fthresh .01)
    endif
    set fmid   = `echo "($fmax+$fthresh)/2.0" | bc -l`
    set fslope = `echo "1.0/($fmax-$fthresh)" | bc -l`
    set cmd = ($cmd -fthresh $fthresh -fmid $fmid -fslope $fslope $maskopt)
    if($Delink) set cmd = ($cmd -delink);
    if($#TCL) set cmd = ($cmd -tcl $TCL)
  else
    # FreeView via tksurferfv
    set cmd = (tksurferfv $trgsubjectuse $hemi $surf -overlay $mapconcatfile)
    if($#annot) set cmd = ($cmd -annot $annot)
    if($surf != white) set cmd = ($cmd -surf white);
    if($surf != pial) set cmd = ($cmd -surf pial);
    if($surf != inflated) set cmd = ($cmd -surf inflated);
    set cmd = ($cmd -fminmax $fthresh $fmax)
  endif

  echo "cd `pwd`"
  echo "setenv SUBJECTS_DIR $SUBJECTS_DIR"
  echo $cmd
  echo ""

  # Print a list of contrasts
  echo "Contrasts"
  @ nth = 0;
  foreach c ($contrastlist)
    echo $nth $c
    @ nth = $nth + 1;
  end
  echo ""
  echo ""

  # Runit
  $cmd

  # delete tmp files
  rm -f $mapconcatfile
  rm -f $sess/$mapfsd/$mapanalysis/$hemi.h.mgh
  rm -f $sess/$mapfsd/$mapanalysis/$hemi.h-offset.mgh

end # Loop over sessions

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

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

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

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

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

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

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

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

    case "-fieldsign":
      set contrastlist = fieldsign;
      set map = fieldsign;
      set DoFieldSign = 1;
      set fthresh = 0.5
      breaksw

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

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

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

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

    case "-no-annot":
    case "-no-aparc":
      set annot = ();
      breaksw

    case "-aparc":
      set annot = aparc
      breaksw

    case "-a2009s":
      set annot = aparc.a2009s
      breaksw

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

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

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

    case "-freeview":
    case "-fv":
      set UseFreeview = 1;
      breaksw

    case "-tks":
    case "-tksurfer":
      set UseFreeview = 0;
      breaksw

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

    case "-tmp":
    case "-tmpdir":
    case "-tmp-dir":
      if($#argv == 0) goto arg1err;
      setenv FREESURFER_TMP_DIR $argv[1];shift;
      breaksw

    case "-nohdr":
      set ShowHDR = 0;
      breaksw

    case "-delink":
      set Delink = 1;
      breaksw

    case "-no-vol2surf":
      set UseVol2Surf = 0;
      breaksw

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

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

    case "-self":
      set trgsubject = self
      breaksw

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

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

    case "-g":
    case "-s":
    case "-sf":
    case "-d":
    case "-df":
      shift;
      # ignore getsesspath arguments 
      breaksw

    case "-cwd":
      # ignore getsesspath arguments 
      breaksw

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

end

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

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

  if($#analysis == 0) then
    echo "ERROR: no analysis specified" 
    exit 1;
  endif
  if(! -e $analysis) then
    echo "ERROR: analysis $analysis does not exist" 
    exit 1;
  endif
  if($#mapanalysis == 0) set mapanalysis = $analysis
  if(! -e $mapanalysis) then
    echo "ERROR: analysis $mapanalysis does not exist" 
    exit 1;
  endif

  ## Parse the analysis
  if($#analysis != 0) then
    if(! -e $analysis ) then
      echo "ERROR: $analysis does not exist." |& tee -a $LF
      exit 1;
    endif
    set infofile = $analysis/analysis.info
    set cfg = $analysis/analysis.cfg
    set fsd = `getana -a $analysis -t fsd`;
    set IsSurface = `getana -a $analysis -t IsSurface`
    if($IsSurface) then
      set RawSpaceType = surface
      set srcsubject = `getana -a $analysis -t subject`
      set RawSpace = $srcsubject;
      set hemi    = `getana -a $analysis -t hemi`
      echo "Surface $srcsubject $hemi"
    endif
    set UseTalairach = `getana -a $analysis -t IsMNI305`;
    if($UseTalairach) then
      set RawSpaceType = volume
      set srcsubject = fsaverage
      set RawSpace = mni305
    endif
    set IsFIR = `getana -a $analysis -t IsFIR`;

    set IsRet = `getana -a $analysis -t IsRetinotopy`
    if($DoFieldSign && ! $IsRet) then
      echo "ERROR: $analysis is not retinotopy, cannot use -fieldsign"
      exit 1;
    endif
    if($IsRet && ! $DoFieldSign && $#contrastlist == 0) then
      set contrastlist = (eccen polar)
      if($#map == 0) set map = fsig;
    endif

    set IsABB = `getana -a $analysis -t IsABBlocked`
    if($IsABB && $#contrastlist == 0) then
      set contrastlist = (fund);
      if($#map == 0) set map = fsig;
    endif
  endif
  if(! $IsFIR) set ShowHDR = 0;
  if($#hemi == 0) then
    echo "ERROR: no hemisphere specified" 
    exit 1;
  endif
  if($#map == 0) set map = sig;

  ## Parse the map analysis
  if($#mapanalysis != 0) then
    set mapfsd = `getana -a $mapanalysis -t fsd`;
    set IsSurface = `getana -a $mapanalysis -t IsSurface`
    if($IsSurface) then
      set MapRawSpaceType = surface
      set mapsrcsubject = `getana -a $mapanalysis -t subject`
      set mapRawSpace = $srcsubject;
      set maphemi    = `getana -a $mapanalysis -t hemi`
      echo "Map Surface $mapsrcsubject $maphemi"
      if($srcsubject != $mapsrcsubject || $hemi != $maphemi) then
        echo "ERROR: map analysis $mapanalysis uses a different surface."
        exit 1;
      endif      
    endif
    set MapUseTalairach = `getana -a $mapanalysis -t IsMNI305`;
    if($UseTalairach != $MapUseTalairach) then
      if($MapUseTalairach) then
        echo "ERROR: map analysis $mapanalysis is mni305-based."
        echo "This does not match the analysis $analysis"
      else
        echo "ERROR: map analysis $mapanalysis is NOT mni305-based."
        echo "This does not match the analysis $analysis"
      endif      
      exit 1;
    endif
  endif

  if($DoAllContrasts) set contrastlist = `getana -a $mapanalysis -t contrasts`
  if($#contrastlist == 0) then
    echo "ERROR: no contrast specified" 
    exit 1;
  endif

  if($#regfile == 0) set regfile = register.dof$RegDOF.dat

  if($ShowHDR && $UseFreeview) then
    echo "ERROR: you must use tksurfer to display the hemodynamic response time course"
    echo "  Two things you can do:"
    echo "    1. To display with freeview but no HRF, re-run with the -nohdr option"
    echo "    2. To display with tksurfer and the HRF, re-run with the -tksurfer option"
  endif

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

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

############--------------##################
usage_exit:
  echo "tksurfer-sess"
  echo ""
  echo "   -s sessid  : session "
  echo "   -hemi hemi : lh or rh";
  echo "   -analysis  analysisname : name of session-level functional analysis";
  echo "   -contrast  contrastname : contrast name"
  echo "   -call      : do all contrasts"
  echo "   -map       mapname      : <sig>, fsig, minsig, iminsig, t"
  echo "   -mapanalysis analysisname : name of session-level functional analysis";
  echo "   -surf surface : default is inflated"
  echo ""
  echo "   -no-aparc : do not show automatic parcellation"
  echo "   -a2009s : use Destrieux atlas instead of Desikian/Killiany"
  echo "   -fthresh fthresh : threshold for color overlay"
  echo "   -fmax    fmax    : saturation threshold for color overlay"
  echo "   -interp method : interpolation for volume analyses"
  echo "   -fwhm fwhm : surface smooth by fwhm (only necessary for volume analyses)"
  echo "   -tcl script.tcl  : run tksurfer with tcl script"
  echo "   -freeview : display with freeview (cannot display HRFs)"
  echo "   -tksurfer  : display with tksurfer"
  echo ""
  echo "   -d  sessdir  ..."
  echo ""

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

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