#! /bin/csh -f

#
# avgraw
#
# average raw fMRI on a time-point by time-point basis
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: nicks $
#    $Date: 2007/01/09 22:41:16 $
#    $Revision: 1.2 $
#
# 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: avgraw,v 1.2 2007/01/09 22:41:16 nicks Exp $'

set n = `echo $argv | grep version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif

set instems  = ();
set outstem  = ();
set stdstem  = ();
set zstem  = ();
set firstslice =  ();
set nslices    = ();
set BFileOutExt = bfloat;
set BFileInExt  = ();
set monly = 0;
set SynthData = 0;
set SynthSeed = -1;

if ( $#argv == 0 ) then
  echo "USAGE: avgraw [-options] -i instem -i instem -o outstem";
  echo "   instem   - prefix of input files";
  echo "Options:";
  echo "   -std stdstem  : save std dev to stdstem"
  echo "   -z   zstem    : save z-score to zstem"
  echo "   -bfloat       : produce bfloat files instead of bshort"
  echo "   -monly mfile  : don't run, just generate a matlab file"
  echo "   -version      : print version and exit"
  echo ""
  echo "  $VERSION"
  echo "  Comments or questions: analysis-bugs@nmr.mgh.harvard.edu"
  echo ""
  exit 1;
endif

goto parse_args;
parse_args_return:

## Autodetect Input extension ##
if($#BFileInExt == 0) then
  set BFileInExt = `getbext $instems[1]`;
  if($status) then
    echo "ERROR: detecting extension of $instems[1]"
    exit 1;
  endif
endif

## Autodetect First Slice ##
if($#firstslice == 0) then
  set firstslice = `getfirstsliceno $instems[1]`;
  if($status) then
    echo "ERROR: detecting slice $firstslice of $instems[1]"
    exit 1;
  endif
endif

## Autodetect Number of Slices ##
if($#nslices == 0) then
  set nslices = `getnslices $instems[1] $firstslice`;
  if($status) then
    echo "ERROR: nslices of $instems[1]"
    exit 1;
  endif
endif

goto check_params;
check_params_return:

set MATLAB = `getmatlab`;
if($status) exit 1;

#### Output Directory #######
if($#outstem > 1) then
  set OutDir  = `dirname  $outstem`;
  set OutBase = `basename $outstem`;
  mkdir -p $OutDir;
endif

if ($monly) then
  set QuitOnError = 0;
  set TARGET = "tee $mfile"
  rm -f $mfile;
else
  set QuitOnError = 1;
  set TARGET =  "$MATLAB -display iconic"
endif  

#---------------------------------------------------------------#
$TARGET <<EOF

  global QuitOnError;
  QuitOnError = $QuitOnError;
  InStems     = splitstring('$instems');
  inext       = '$BFileInExt';
  outstem     = '$outstem';
  stdstem     = '$stdstem';
  zstem       = '$zstem';
  outext      = '$BFileOutExt';
  FirstSlice  = $firstslice;
  NSlices     = $nslices;
  SynthData   = $SynthData;
  SynthSeed   = $SynthSeed;

  if(SynthSeed == -1) SynthSeed = sum(100*clock); end
  fprintf('SynthSeed = %g\n',SynthSeed);
  randn('state',SynthSeed); 

  LastSlice = FirstSlice + NSlices - 1;
  NInputs = size(InStems,1);

  for slice = FirstSlice:LastSlice
    ysum  = 0;
    ysum2 = 0;
    for input = 1:NInputs
      instem = deblank(InStems(input,:));
      infile = sprintf('%s_%03d.%s',instem,slice,inext); 
      fprintf('Loading %s\n',infile);
      y = fmri_ldbfile(infile);
      if($SynthData) y = randn(size(y)); end
      ysum = ysum + y;
      if(~isempty(stdstem) | ~isempty(zstem))
        ysum2 = ysum2 + y.^2; 
      end
    end
    yavg = ysum/NInputs;
    outfile = sprintf('%s_%03d.%s',outstem,slice,outext); 
    fprintf('Saving %s\n',outfile);
    fmri_svbfile(yavg,outfile);

    if(~isempty(stdstem) | ~isempty(zstem))
      ystd = sqrt(ysum2/(NInputs-1) - yavg.^2);
    end

    if(~isempty(stdstem))
      outfile = sprintf('%s_%03d.bfloat',stdstem,slice);
      fprintf('Saving %s\n',outfile);
      fmri_svbfile(ystd,outfile);
    end

    if(~isempty(zstem))
      izero = find(ystd == 0);
      ystd(izero) = 10^5;
      z = yavg./ystd;
      outfile = sprintf('%s_%03d.bfloat',zstem,slice);
      fprintf('Saving %s\n',outfile);
      fmri_svbfile(z,outfile);
    end

  end
  return;

  if(QuitOnError)  quit; end
  
EOF

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 instems = ($instems $argv[1]); shift;
      breaksw

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

    case "-o":
      if ( $#argv == 0) goto arg1err;
      if ( $#outstem != 0 ) then
        echo ERROR: only one outstem allowed.
        exit 1
      endif
      set outstem = $argv[1]; shift;
      breaksw

    case "-std":
      if ( $#argv == 0) goto arg1err;
      if ( $#stdstem != 0 ) then
        echo ERROR: only one stdstem allowed.
        exit 1
      endif
      set stdstem = $argv[1]; shift;
      breaksw

    case "-z":
      if ( $#argv == 0) goto arg1err;
      if ( $#zstem != 0 ) then
        echo ERROR: only one zstem allowed.
        exit 1
      endif
      set zstem = $argv[1]; shift;
      breaksw

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

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

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

    case "-bfloat":
      set BFileOutExt  = bfloat;
      breaksw

    case "-bshort":
      set BFileOutExt  = bfloat;
      breaksw

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

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

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

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

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

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

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

############--------------##################
check_params:
  if ($#instems < 2) then
     echo "ERROR: must specify at least 2 input stems";
     exit 1
  endif

  if ($#outstem == 0) then
     echo "ERROR: must specify an output stem";
     exit 1
  endif

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


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