#! /bin/csh -f

#
# stc.fsl - front-end for FSL's slicetimer slice timing correction
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: greve $
#    $Date: 2011/03/28 20:38:31 $
#    $Revision: 1.7.4.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: stc.fsl,v 1.7.4.1 2011/03/28 20:38:31 greve Exp $'
set debug = 0;
set cmdargs = ($argv);
set PrintHelp = 0;

set infunc = ();
set outfunc = ();
set SliceSeq = ();
set OutDataType = ();
set OutScale = ();

set tmpdir = ();
set cleanup = 1;


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

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:


set outdir = `dirname $outfunc`;
if($#tmpdir == 0)  set tmpdir = $outdir/tmp.$$

if($tmpdir == $outdir && $cleanup) then
  echo "ERROR: tmpdir and outdir are the same ($tmpdir)"
  echo "  If you really want to do this, run with --nocleanup"
  exit 1;
endif

mkdir -p $outdir
mkdir -p $tmpdir

set LF = $outfunc.log
if(-e $LF) mv $LF $LF.bak

date | tee -a $LF
pwd  | tee -a $LF
echo $0  | tee -a $LF
echo $cmdargs | tee -a $LF
echo $VERSION | tee -a $LF
hostname  | tee -a $LF

# Get number of input slices
set tmpfile = $tmpdir/nslices.dat
set cmd = (mri_info --nslices --o $tmpfile $infunc);
pwd |& tee -a $LF
echo $cmd  |& tee -a $LF
$cmd  |& tee -a $LF
if($status) exit 1;
set nslices = `cat $tmpfile`;
echo "nslices =  $nslices" |& tee -a $LF

# Choose the proper slice sequence for siemens default
if($SliceSeq == siemens) then
  # Determine whether the number of slices is odd
  set isodd = `echo "$nslices%2" | bc`;
  echo "isodd = $isodd" |& tee -a $LF
  if($isodd)   set SliceSeq = odd
  if(! $isodd) set SliceSeq = even
  echo "Changing siemens slice seq to $SliceSeq"|& tee -a $LF
endif

# slice acq order file
set saof = $tmpdir/saof.dat 
rm -f $saof

# Construct the slice sequence
switch($SliceSeq)
  case "odd":
    @ nthslice = 1
    while($nthslice <= $nslices)
      echo $nthslice >> $saof
      @ nthslice = $nthslice + 2;
    end
    @ nthslice = 2
    while($nthslice <= $nslices)
      echo $nthslice >> $saof
      @ nthslice = $nthslice + 2;
    end
    breaksw
  case "even":
    @ nthslice = 2
    while($nthslice <= $nslices)
      echo $nthslice >> $saof
      @ nthslice = $nthslice + 2;
    end
    @ nthslice = 1
    while($nthslice <= $nslices)
      echo $nthslice >> $saof
      @ nthslice = $nthslice + 2;
    end
    breaksw
  case "up":
    @ nthslice = 1
    while($nthslice <= $nslices)
      echo $nthslice >> $saof
      @ nthslice = $nthslice + 1;
    end
    breaksw
  case "down":
    @ nthslice = $nslices
    while($nthslice > 0)
      echo $nthslice >> $saof
      @ nthslice = $nthslice - 1;
    end
    breaksw
endsw

echo "---------- Slice Sequence ---------------" | tee -a $LF
cat $saof | tee -a $LF
echo "-----------------------------------------" | tee -a $LF

# Convert the input to analyze
set intmp = $tmpdir/in.nii
set cmd = (mri_convert $infunc $intmp)
pwd |& tee -a $LF
echo $cmd  |& tee -a $LF
$cmd  |& tee -a $LF
if($status) exit 1;
set intmpmat = $tmpdir/in.mat

# Apply slicetimer
setenv FSLOUTPUTTYPE NIFTI
set outtmp = $tmpdir/out.nii
set cmd = (slicetimer.fsl -i $intmp -o $outtmp --ocustom=$saof)
pwd |& tee -a $LF
echo $cmd  |& tee -a $LF
$cmd  |& tee -a $LF
if($status) exit 1;

if(! -e $outtmp) then
  echo "ERROR: failed to create $outtmp" |& tee -a $LF
  exit 1;
endif
set outtmpmat = $tmpdir/out.mat

# Copy the matfile
# cp $intmpmat $outtmpmat 

# Convert back to the output
set cmd = (mri_convert $outtmp $outfunc)
if($#OutDataType) set cmd = ($cmd -odt $OutDataType)
if($#OutScale)    set cmd = ($cmd --scale $OutScale)

pwd |& tee -a $LF
echo $cmd  |& tee -a $LF
$cmd  |& tee -a $LF
if($status) exit 1;

if($cleanup) rm -rf $tmpdir

echo "stc.fsl Done"    | tee -a $LF

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


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

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

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

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

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

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

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

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

    case "--odd":
      set SliceSeq = "odd"; 
      breaksw
    case "--even":
      set SliceSeq = "even"; 
      breaksw
    case "--up":
      set SliceSeq = "up"; 
      breaksw
    case "--down":
      set SliceSeq = "down"; 
      breaksw
    case "--siemens":
      set SliceSeq = "siemens"; 
      breaksw

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

end

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

############--------------##################
check_params:
if($#infunc == 0) then
  echo "ERROR: need input functional"
  exit 1;
endif
if($#outfunc == 0) then
  echo "ERROR: need output functional"
  exit 1;
endif
if(! -e $infunc) then
  echo "ERROR: input functional $infunc does not exist"
  exit 1;
endif
if($#SliceSeq == 0) then
  echo "ERROR: must supply a slice sequence "
  exit 1;
endif
goto check_params_return;
############--------------##################

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

############--------------##################
usage_exit:
  echo ""
  echo "USAGE: stc.fsl"
  echo ""
  echo "   --i infunc "
  echo "   --o outfunc "
  echo "   --up   : sequential starting on first slice"
  echo "   --down : sequential starting on last slice"
  echo "   --odd  : interleaved starting on first slice"
  echo "   --even : interleaved starting on second slice"
  echo "   --siemens : use siemens default slice timing" 

  echo "   --debug"
  echo "   --tmp tmpdir     : use tmpdir as temporary directory"
  echo "   --version        : print version and exit"
  echo "   --help           : print help and exit"
  echo ""

  if(! $PrintHelp) exit 1;

  echo $VERSION

  echo "------------------------------------------------------------"
  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

Front-end for FSLs slice timing correction program (slicetimer). The slice
order is the list of slices in the order in which they were acquired. Eg,
a slice order of 1 3 2 4 means that the first slice in the stack stored 
on disk was the first acquired, the third slice in the stack stored 
on disk was the second acquired, etc.

ARGUMENTS

--i input functional

Any input format supported by mri_convert

--o output functional

Any output format supported by mri_convert except analyze3d/spm

--up

Sequential starting on first slice: 1 2 3 4 5 ...

--down

Sequential starting on last slice: nslices (nslices-1) (nslices-2) ... 1

--odd 

Interleaved starting on first slice: 1 3 5 7 ... 2 4 6 8 ...

--even

Interleaved starting on second slice: 2 4 6 8 ... 1 3 5 7 ... 

--siemens

If the number of slices is odd, then uses --odd. If the number of slices is even, 
then uses --even.

NOTES

1. Slice order starts at 1.

2. Delay=0 is set at the "relative center" and is at: (nslices/2 -
   0.5)/nslices. Eg, if nslices=20, then the relative center is .475
   through the TR (which corresponds to half way through the 9th
   acquired slice). If nslices=21, then delay=0 at the 10th slice.
   As a result, you must shift the timing in your stimulus paradigms
   by subtracting half a TR from the onset times.