#! /bin/tcsh -f

#
# reg-feat2anat
#
# Registers example func to freesurfer using fsl
#
# 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 = 'feat2surf dev';
set FeatDirList = ();
set ProjFrac = ();
set HemiList = (lh rh);
set Surf = white;
set TargetSubj = fsaverage;
set debug = 0;
set DontRun = 0;
set usedev = 0;
set myoutstatsdir = ();
set NoLog = 0;
set fslext = ();
set CopeOnly = 0;

set PrintHelp = 0;

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

source $FREESURFER_HOME/sources.csh

# Parse the command-line arguments
goto parse_args;
parse_args_return:

# Check the command-line arguments
goto check_params;
check_params_return:

if($?DEV == 0) set DEV = "";

foreach FeatDir ($FeatDirList)

# Set up a log file
if(! $NoLog) then
  set LF = $FeatDir/reg/freesurfer/feat2surf.log
  echo log file is $LF
  if(-e $LF) mv $LF $LF.bak
else
  set LF = /dev/null
endif
echo FeatDir is $FeatDir | tee -a $LF
date | tee -a $LF
pwd | tee -a $LF
echo $0 | tee -a $LF
echo $cmdargs | tee -a $LF
uname -a | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "subject is $subject" | tee -a $LF

set statsdir = $FeatDir/stats
set reg_anat2exf = $FeatDir/reg/freesurfer/anat2exf.register.dat
if(! -e $reg_anat2exf) then
  echo "ERROR: cannot find $reg_anat2exf"
  echo "Try running reg-feat2anat"
  exit 1
endif
#set subject = `cat $reg_anat2exf| head -n 1`;
set subject = `reg2subject --r $reg_anat2exf`;


# Automatically detect the file format
set exfbase = $FeatDir/example_func
if(-e $exfbase.nii.gz) set fslext = nii.gz;
if(-e $exfbase.nii)    set fslext = nii;
if(-e $exfbase.img)    set fslext = img;
if($#fslext == 0) then
  echo "ERROR: cannot determine type of fsl output"
  exit 1;
endif

# Get a list of images to convert
#set nonomatch # This prevents crashing due to file expansion errors
if($CopeOnly) then
  set imglist = (`ls $statsdir/cope*.$fslext $statsdir/varcope*.$fslext`)
  set subjectlist = ($TargetSubj)
else
  set imglist = (`ls $statsdir/*.$fslext`)
  set imglist = ($imglist `ls $imglist $FeatDir/cluster*.$fslext`)
  set imglist = ($imglist `ls $imglist $FeatDir/rendered_thresh*.$fslext`)
  set subjectlist = ($subject $TargetSubj)
endif

set nimg = $#imglist
if($nimg == 0) then
  echo "ERROR: cannot find any .img, .nii, .nii.gz files in $statsdir"| tee -a $LF
  exit 1;
endif
echo "Found $nimg images in $statsdir" | tee -a $LF

@ nthimg = 0
foreach img ($imglist)
   @ nthimg = $nthimg + 1;
   echo $img ---------------------

   # Get the extension
   set ext = $fslext
   set baseimg = `basename $img $ext`

   # Dont know whether these should be conv too
   #if($baseimg == corrections) continue;
   #if($baseimg == threshac1) continue;

   foreach targ ($subjectlist)
     foreach hemi ($HemiList)

       if($targ != $subject) then
         set outext = $fslext
         set reshape = (--reshape)
       else
         set outext = mgh
         set reshape = ()
       endif

       # Only used for testing
       if($#myoutstatsdir == 0) then
         set outstatsdir = $FeatDir/reg_surf-$hemi-$targ/stats
       else
         set outstatsdir = $myoutstatsdir
       endif

       mkdir -p $outstatsdir
       set outimg = $outstatsdir/$baseimg$outext

       set V2S = mri_vol2surf
       if($usedev) set V2S = $DEV/$V2S/$V2S

      set cmd = ($V2S --src $img --srcreg $reg_anat2exf \
                  --trgsubject $targ --hemi $hemi --out $outimg\
                  --surf $Surf $reshape)
       if($#ProjFrac) set cmd = ($cmd --projfrac $ProjFrac)

       echo " " | tee -a $LF
       echo "#@# $baseimg $targ $hemi $nthimg/$nimg ----=====-----=====----" | tee -a $LF
       date      | tee -a $LF
       pwd       | tee -a $LF
       echo $cmd | tee -a $LF
       if(! $DontRun) $cmd  | tee -a $LF
       if($status) then
         echo "ERROR: with $cmd" | tee -a $LF
         pwd
         exit 1;
       endif
       echo " " | tee -a $LF
     end # hemi list
  end # target subject

  @ nthimg = $nthimg + 1;
end # img list

end # FeatDir List

date | tee -a $LF
echo "feat2surf done" | tee -a $LF

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


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

  set flag = $argv[1]; shift;

  switch($flag)

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

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

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

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

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

    case "--hemi":
      if ( $#argv == 0) goto arg1err;
      set HemiList = $argv[1]; shift;
      if($HemiList != lh && $HemiList != rh) then
        echo "--hemi = $HemiList, must be either lh or rh"
        exit 1;
      endif
      breaksw

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

    case "--cope-only":
      set CopeOnly = 1;
      breaksw

    case "--debug":
      set verbose = 1;
      set echo = 1; # turns on terminal echoing
      set debug = 1;
      breaksw

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

    case "--dontrun":
      set DontRun = 1;
      breaksw

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

    default:
      echo "ERROR: flag $flag not recognized"
      exit 1;
      breaksw
  endsw

end

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

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

if($#FeatDirList == 0) then
  echo "ERROR: must specify at least one feat dir"
  exit 1;
endif

foreach FeatDir ($FeatDirList)
  set reg_anat2exf = $FeatDir/reg/freesurfer/anat2exf.register.dat
  if(! -e $reg_anat2exf) then
    echo "ERROR: cannot find $reg_anat2exf. You must run reg-feat2anat first."
    exit 1;
  endif
  if($?SUBJECTS_DIR == 0) then
    echo "ERROR: FreeSurfer environment variable SUBJECTS_DIR not defined"
    exit 1;
  endif
  if(! -e $SUBJECTS_DIR) then
    echo "ERROR: FreeSurfer SUBJECTS_DIR $SUBJECTS_DIR not found"
    exit 1;
  endif

  set subject = `cat $reg_anat2exf| head -n 1`;
  if(! -e $SUBJECTS_DIR/$subject) then
    echo "ERROR: cannot find FreeSurfer subject $subject in $SUBJECTS_DIR"
    exit 1;
  endif

  foreach hemi ($HemiList)
    set surf = $SUBJECTS_DIR/$subject/surf/$hemi.$Surf
    if(! -e $surf) then
      echo "ERROR: cannot find surface $hemi.$Surf in FreeSurfer $SUBJECTS_DIR/$subject/surf"
      exit 1;
    endif
    if($#ProjFrac) then
      set thickness = $SUBJECTS_DIR/$subject/surf/$hemi.thickness
      if(! -e $thickness) then
        echo "ERROR: cannot find surface $thickness in FreeSurfer $SUBJECTS_DIR/$subject/surf"
        exit 1;
      endif
    endif
  end
end # FeatDirList

#if($?FSLOUTPUTTYPE == 0) set FSLOUTPUTTYPE = ANALYZE
#if($FSLOUTPUTTYPE != ANALYZE && $FSLOUTPUTTYPE != NIFTI) then
#  echo "ERROR: FSLOUTPUTTYPE = $FSLOUTPUTTYPE, must be ANALYZE or NIFTI"
#  exit 1;
#endif


endif

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

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

############--------------##################
usage_exit:
  echo ""
  echo "USAGE: feat2surf"
  echo ""
  echo " --feat dir  <--feat dir>  : Feat output directory"
  echo " --featdirfile file : file with a list of feat directories"
  echo ""
  echo "Optional flags and arguments:"
  echo ""
  echo " --projfrac fraction : project fraction into cortical surface (default is 0)"
  echo " --hemi hemi : run on hemi (lh or rh) only. Default is lh and rh. "
  echo " --target subject : subject used to define common surface space (default is fsaverage)"
  echo " --surf surface : surface to resample to (default is white)"
  echo " --cope-only : only map the copes and varcopes to the common surface space"
  echo " "
  echo " --version : print version and exit"
  echo " --help    : print help and exit"
  echo " --debug   : turn on debugging"
  echo " --nolog   : do not create a log file"
  echo " --out outdir : put output in outdir instead of feat/reg_surf-?h/stats (debugging)"
  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

Resamples Feat statistics onto the surface of the subject and onto a
stereo-taxic surface atlas.

Statistics are obtained from FeatDir/stats as well as those volumes in
FeatDir beginning with "cluster" and "rendered". The results are
stored in FeatDir/reg_surf-lh-SUBJECT/stats, and
FeatDir/reg_surf-rh-SUBJECT/stats where SUBJECT is the freesurfer name
of the subject. It will also resample to the common space and stored
the results in in FeatDir/reg_surf-lh-targid/stats and
FeatDir/reg_surf-rh-targid/stats, where targid is the subject that
defines the common surface space (default is fsaverage). If you are using
the "new" FreeSurfer recon where you have registered your surfaces
to the "filled" atlas, then use fsaverage instead.

Works properly on ANALYZE or NIFTI; output is stored in mgz format.

Creates log file featdir/reg/freesurfer/feat2surf.log

REQUIRED ARGUMENTS

--feat featdir <--feat featdir ...>

Directory where Feat results are stored. Can specify multiple --feat.

--featdirfile file <--featdirfile file>

File with a list of Feat directories. Both --feat and --featdirfile
can be specified (multiple times).

OPTIONAL ARGUMENTS

--projfrac fraction

Sample functional a fraction of the cortical thickness normal to the
surface. Eg, --projfrac 0.5 will sample in the middle of the cortical
ribbon. Default is 0.  The subject must have thickness files in the
surf directory.

--hemi hemi

Run on hemi (lh or rh) only. Default is to do both.

--target subject

Subject used to define common surface space (default is fsaverage).

--surf surface

Surface to resample to (default is white).

EXAMPLE

# Assume a data set exists called fbert.img
# Analyze fbert.img with Feat to generate directory fbert.feat

# Register fbert.feat to subject bert:
reg-feat2anat --feat fbert.feat --subject bert

# Resample fbert.feat onto the surface
feat2surf --feat fbert.feat

# View on native surface
tksurfer bert lh inflated \
  -overlay fbert.feat/reg_surf-lh-bert/stats/zstat1.nii.gz

# View on common surface
tksurfer fsaverage lh inflated \
  -overlay fbert.feat/reg_surf-lh-fsaverage/stats/zstat1.nii.gz


BUGS

The surface-based data are stored in "volume" formats (ie, the format
found in the FEAT dir).  In this case, the spatial dimensions do not
correspond to volume dimensions but rather the number of vertices on
the surface. This is typically on the order of 150,000. Ideally, we
would simply store the surface data as a "volume" with 150,000
columns, 1 row, and 1 slice.  However, ANALYZE cannot represent
dimensions over 2^15, or 32768.  Same problem with NIFTI. To get
around this, we factor the number of vertices into dimensions smaller
than 32768. Eg, the fsaverage surface has 163842 vertices which we
factor into 1974x83. However, for individuals there is no guarantee
that the number of vertices will have a factor less than 32768. For
this reason the data resampled on the individual surfaces is stored in
MGH format (which does not have these restrictions).