#! /bin/csh -f
#
# fsfeat - freesurfer front end for FSL's FEAT
# For more info on FEAT see www.fmrib.ox.ac.uk/fsl

# Things to do:
#   2. init mtxs
#   5. compare TR


#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: greve $
#    $Date: 2007/06/27 06:13:53 $
#    $Revision: 1.14 $
#
# 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: fsfeat,v 1.14 2007/06/27 06:13:53 greve Exp $';
set inputargs = ($argv);

set fsf     = ();
set invol   = ();
set schdir  = ();
set outdir  = ();
set subject = ();
set keep_residuals = 0;
set UsePatch = 0;

set PrintHelp = 0;

#### If no arguments, print usage and exit ####
if($#argv == 0) then
  goto usage_exit;
  exit 1;
endif
set n = `echo $argv | grep -e -help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
  exit 1;
endif

##### Print out version info, if needed ####
set n = `echo $argv | grep version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

if($FSLOUTPUTTYPE == NIFTI)    set outext = nii;
if($FSLOUTPUTTYPE == NIFTI_GZ) set outext = nii.gz
if($FSLOUTPUTTYPE == ANALYZE) then
  echo "ERROR: FSLOUTPUTTYPE == ANALYZE, get a real format!"
  exit 1;
endif

# Set up the output path
#   Feat needs a full path
#   But this cannot exist, or else get +s
set outdirpath = $outdir.feat
mkdir -p $outdirpath
set outdirpath = `getfullpath $outdirpath`;
rm -fr $outdirpath
echo $outdir

# Get number of time points from input volume
avwinfo $invol | awk '{if($1 == "dim4") print $2}'
if($status) exit 1;
set nframes = `avwinfo $invol |& awk '{if($1 == "dim4") print $2}'`

# Create initial example_func to std matrix
set stdvol = $FSLDIR/etc/standard/avg152T1.img
set initmtx = /tmp/fsfeat.design.$$.initfsl.mtx
set cmd = (tkregister2_cmdl --targ $stdvol --mov $invol \
   --reg /tmp/$$.reg --fslregout $initmtx --regheader --noedit);
echo $cmd 
$cmd
if($status) exit 1;

# Temporary fsf file
set tmpfsf = /tmp/fsfeat.design.$$.fsf
cp $fsf $tmpfsf
echo $tmpfsf

# Add customizations
echo "#-------------------------------------------" >> $tmpfsf
echo "#  fsfeat customizations"                     >> $tmpfsf
echo "#-------------------------------------------" >> $tmpfsf
echo "# `date`"   >> $tmpfsf
echo "# $VERSION" >> $tmpfsf
echo "set feat_files(1)" '"'$invol'"' >> $tmpfsf
echo "set fmri(init_standard)" '"'$initmtx'"' >> $tmpfsf
echo "set fmri(npts) $nframes" >> $tmpfsf
# always keep residuals here; maybe deleted later
echo "set fmri(cleanup_residuals_yn) 0" >> $tmpfsf
# init_initial_highres init_initial_highres

set nEVs = `fsfget $fsf evs_orig`;
@ nthEV = 0;
while ($nthEV < $nEVs)
  @ nthEV = $nthEV + 1;
  set tmp = `grep fmri\(custom$nthEV\) $fsf | awk '{print $3}'| sed 's/\"//g'`
  if($#tmp == 0) continue;
  if($#schdir == 0)  set schdir = `dirname $invol`;
  set sch = `basename $tmp`;
  echo $nthEV $sch
  set schpath = $schdir/$sch
  if(! -e $schpath) then
    echo "ERROR: cannot find $schpath"
    exit 1;
  endif
  echo "set fmri(custom$nthEV)" '"'$schpath'"' >> $tmpfsf
end

echo "set fmri(outputdir)" '"'$outdirpath'"'  >> $tmpfsf
echo "set fmri(featwatcher_yn) 0"             >> $tmpfsf

if(! $UsePatch) feat $tmpfsf
if($UsePatch)   feat.fsfast-patch $tmpfsf
if($status) exit 1;

set exf = $outdirpath/example_func.$outext
if(! -e $exf) then
  echo "ERROR: feat failed"
  exit 1;
endif

set mask = $outdirpath/mask

# Compute the mean image (after intensity scaling)
set mean_func = $outdirpath/mean_func
set cmd = (avwmaths_32R $outdirpath/filtered_func_data \
  -Tmean $mean_func)
pwd
echo $cmd
$cmd
if($status) exit 1;

# Compute residual variance
set res4d = $outdirpath/stats/res4d
set rvar  = $outdirpath/stats/rvar
set dof = `cat $outdirpath/stats/dof`
set cmd = (avwmaths_32R  $res4d -sqr -Tmean -div $dof -mul $nframes $rvar)
pwd
echo $cmd
$cmd
if($status) exit 1;

# Compute residual stddev
set rstd  = $outdirpath/stats/rstd
set cmd = (avwmaths_32R  $rvar -sqrt $rstd)
pwd
echo $cmd
$cmd
if($status) exit 1;

# Compute FSNR
set fsnr  = $outdirpath/stats/fsnr
set cmd = (avwmaths_32R  $mean_func -div $rstd $fsnr)
pwd
echo $cmd
$cmd
if($status) exit 1;

# Compute FSNR2
set fsnr2  = $outdirpath/stats/fsnr2
set cmd = (avwmaths_32R  $fsnr -sqr $fsnr2)
pwd
echo $cmd
$cmd
if($status) exit 1;

# Compute averages over the mask
set qadir = $outdirpath/qa
mkdir -p $qadir
avwmeants -m $mask -i $exf       -o $qadir/example_func.dat
avwmeants -m $mask -i $mean_func -o $qadir/mean_func.dat # always 1000
avwmeants -m $mask -i $rvar      -o $qadir/rvar.dat
avwmeants -m $mask -i $rstd      -o $qadir/rstd.dat
avwmeants -m $mask -i $fsnr      -o $qadir/fsnr.dat
avwmeants -m $mask -i $fsnr2     -o $qadir/fsnr2.dat

# Number of voxels in the mask
cat $outdirpath/stats/smoothness | \
    awk '{if($1 == "VOLUME") print $2}' > $qadir/nmask.dat

# Redo smoothness estimation, saving the fwhm
set cmd = (smoothest -d $dof -r $res4d -m $mask -V)
echo $cmd
$cmd >& $qadir/fwhm.log
if($status) exit 1;

set fwhmx = `grep FWHM $qadir/fwhm.log | grep mm | awk '{print $3}'`;
set fwhmy = `grep FWHM $qadir/fwhm.log | grep mm | awk '{print $7}'`;
set fwhmz = `grep FWHM $qadir/fwhm.log | grep mm | awk '{print $11}'`;
echo "$fwhmx $fwhmy $fwhmz" > $qadir/fwhm.dat

# Now put it all in one file
rm -f $qadir/qa.dat
foreach q (nmask mean_func rstd fsnr fwhm example_func fsnr2 rvar )
  set v = `cat $qadir/$q.dat`;
  echo -n "$v ">> $qadir/qa.dat
end
echo "" >> $qadir/qa.dat

# Now delete the residuals
if(! $keep_residuals) rm $res4d

# Compute sigs from zstats
set statsdir = $outdirpath/stats
set nzstats = `ls $statsdir/zstat*.$outext | wc -w`;
echo nzstats = $nzstats
set concatlist = ()
@ nthzstat = 0;
while ($nthzstat < $nzstats)
  @ nthzstat = $nthzstat + 1;
  set z = $statsdir/zstat$nthzstat.$outext
  set sig = $statsdir/zsig$nthzstat.$outext
  set cmd = (mri_z2p --z $z --log10p $sig --two-sided) #unsigned
  pwd
  echo $cmd
  $cmd
  if($status) exit 1;
  # symlink to meaningful name here
  set conname = `fsfget $fsf conname_real.$nthzstat`
  if($#conname != 0) then
    pushd $statsdir > /dev/null
    set tmp = "$conname.sig.$outext"
    if(-e "$tmp") rm "$tmp";
    ln -s zsig$nthzstat.$outext "$tmp";
    popd  > /dev/null
  endif
  set concatlist = ($concatlist $sig);
end
# concatenate them
mri_concat $concatlist --o $statsdir/zsig.all.$outext

# Compute sigs from zfstats
set nzfstats = `ls $statsdir/zfstat*.$outext | wc -w`;
echo nzfstats = $nzfstats
@ nthzfstat = 0;
while ($nthzfstat < $nzfstats)
  @ nthzfstat = $nthzfstat + 1;
  set zf  = $statsdir/zfstat$nthzfstat.$outext
  set sig = $statsdir/zfsig$nthzfstat.$outext
  set cmd = (mri_z2p --z $zf --log10p $sig --one-sided) # signed
  pwd
  echo $cmd
  $cmd
  if($status) exit 1;
end

if($#subject) then
  set cmd = (reg-feat2anat --feat $outdirpath --subject $subject);
  echo "setenv SUBJECTS_DIR $SUBJECTS_DIR"
  echo "cd `pwd`"
  echo $cmd
  $cmd
  if($status) exit 1;
endif

echo " " 
echo "fsfeat done" 
echo " " 

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

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

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

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

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

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

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

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

    case "--svres":
      set keep_residuals = 1;
      breaksw

    case "--patch":
      set UsePatch = 1;
      breaksw

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

    default:
      echo "ERROR: flag $flag unrecognized"
      breaksw
  endsw

end

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

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

  if($#fsf == 0) then
    echo "ERROR: must specify a design.fsf"
    exit 1;
  endif

  if($#invol == 0) then
    echo "ERROR: must specify a input volume"
    exit 1;
  endif
  if(! -e $invol) then
    echo "ERROR: cannot find $invol"
    exit 1;
  endif
  isanalyze $invol
  if($status) then
    echo "ERROR: get a real format"
    exit 1;
  endif
  set invol = `getfullpath $invol`;

  if(! -e $fsf) then
    echo "ERROR: cannot find $fsf"
    exit 1;
  endif

  if($#schdir) then
    if(! -e $schdir) then
      echo "ERROR: cannot find $schdir"
      exit 1;
    endif
    pushd $schdir > /dev/null
    set schdir = `pwd`;
    popd > /dev/null
  endif

  if($#subject) then
    set sd = $SUBJECTS_DIR/$subject
    if(! -e $sd) then
      echo "ERROR: cannot find $subject in $SUBJECTS_DIR"
      exit 1;
    endif
  endif

  if($#outdir == 0) then
    echo "ERROR: must specify an outdir"
    exit 1;
  endif
  if(-e $outdir) then
    pushd $outdir > /dev/null
    set outdir = `pwd`;
    popd > /dev/null
  endif

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

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

############--------------##################
usage_exit:
  echo ""
  echo "fsfeat "
  echo ""
  echo "  --fsf     design.fsf"
  echo "  --i       input volume"
  echo "  --schdir  schedule directory (if needed)"
  echo ""
  echo "  --outdir  outputdir : will be outputdir.feat"
  echo ""
  echo "  --subject subject (optional)"
  echo "  --svres   : save resduals"
  echo ""
  echo "  --patch : use FSFAST_HOME/bin/feat.patch"
  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

Help!