#! /bin/csh -f

#
# synth-phenc
#
# fMRI phase-encoding synthsizer
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: nicks $
#    $Date: 2007/01/09 22:41:19 $
#    $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: synth-phenc,v 1.2 2007/01/09 22:41:19 nicks Exp $';
set cmdargv = "$argv";

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

set sfastem  = ();
set outstem = ();
set bext = bshort;
set monly = 0;

set TR = 4;
set offset = 1000;
set trend  = 0;
set ncycles = 8;
set direction = "pos";
set nrows =  64;
set ncols =  64;
set ntp   = 128;
set firstslice = 0;
set nslices = 1;
set phase = 0;
set snr = 1;
set pctsigch = 5;
set scale = 1;
#set roi = (20 20 40 40);
set roi = ();
set delay      = 0;
set dispersion = -1;  # default to sinusoid #
set slice_delay_file = ();
set nskip = 0;
set seed = ();

if($#argv == 0) goto usage_exit;

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

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

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

set outdir = `dirname $outstem`;
mkdir -p $outdir

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

  global QuitOnError;
  QuitOnError = $QuitOnError;
  outstem    = '$outstem';
  bext       = '$bext';

  TR         = $TR;
  firstslice = $firstslice;
  nslices    = $nslices;
  offset     = [$offset];
  trend      = [$trend];
  ncycles    = [$ncycles];
  direction  = '$direction';
  nrows      = [$nrows];
  ncols      = [$ncols];
  ntp        = [$ntp];
  ph         = [$phase]*pi/180;
  snr        = [$snr];
  pctsigch   = [$pctsigch];
  roi        = [$roi];
  delay      = [$delay];
  dispersion = [$dispersion];
  scale      = [$scale];
  nskip      = [$nskip];
  seed       = [$seed];

  %% Set the seed of the random number generator %%
  if(isempty(seed)) seed = sum(100*clock); end
  randn('state',seed);
  rand('state',seed);
  fprintf('Seed = %d\n',seed);

  Trun = TR*ntp;
  fund = ncycles/Trun;
  Nv = nrows * ncols;
  % Nv = 1;

  %% Modify the phase for negative direction %%
  if( strcmp(direction,'neg')) 
     ph = -ph; 
  end

  if(~isempty(roi))
    tmp = zeros(nrows,ncols);
    tmp(roi(1):roi(3),roi(2):roi(4)) = 1;
    indroi  = find(tmp==1);
    indnone = find(tmp==0);
  else
    indroi  = [];
    indnone = [];
  end

  lastslice = firstslice + nslices - 1;
  for slice = firstslice:lastslice,

    if(snr ~= 0)
      signal = fmri_synth_phenc_sig(fund, TR, ntp, ph, delay, dispersion);
      psc    = 100*max(signal)./offset;
      signal = signal*pctsigch./psc;
      signalpower = (std(signal)).^2;
      signal = repmat(signal, [1 Nv]);
      if(~isempty(indnone)) signal(:,indnone) = 0; end
    else
      signalpower = 1;
    end

    if(snr < 0)
      % No noise %
      noise = 0;
    elseif(snr == 0)
      % All noise %
      signal = 0;
      noise  = randn(ntp,Nv) ;
      psc    = 100*max(noise,[],1)./offset;
      psc    = repmat(psc, [ntp 1]);
      noise  = noise*pctsigch./psc;
    else % snr > 0 %
      % Noise and Signal %
      % this will generate an error %
      noise = sqrt(signalpower) * randn(ntp,Nv)/snr  ;
    end

    if(trend ~= 0)
      t = [0:ntp-1]'; %'
      trendsig = repmat(trend*t, [1 Nv]);
    else
      trendsig = 0;
    end

    data = scale * (signal + noise + offset + trendsig);

    if(nskip > 0) 
      data(1:nskip,:) = 0;
    end

    data = reshape(data', [nrows ncols ntp]); %'
   
    fname = sprintf('%s_%03d.%s',outstem,slice,bext);
    fprintf('Saving to %s\n',fname);
    fmri_svbfile(data,fname);

  end

  if(QuitOnError)  quit; end
  
EOF

exit 0;

############--------------##################
parse_args:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    case "-snr":
      if ( $#argv == 0) goto arg1err;
      set snr = $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 "-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 "-ntp":
      if ( $#argv == 0) goto arg1err;
      set ntp = $argv[1]; shift;
      breaksw

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  if($direction != "pos" && $direction != "neg") then
    echo "ERROR: direction must be either 'pos' or 'neg'";
    exit 1;
  endif    

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


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


############--------------##################
usage_exit:
  echo "USAGE: synth-phenc -o outstem"
  echo "  -TR TR  : temporal resolution of synthesized scan (sec)"
  echo "  -offset offset   : offset at each voxel (MR units)"
  echo "  -trend  trend    : slope of temporal trend"
  echo "  -ncycles ncycles : number of stimulation cycles"
  echo "  -nskip   nskip   : set first nskip images to zero"
  echo "  -direction dir   : <pos> or neg"
  echo "  -phase     ph    : degrees"
  echo "  -delay     del   : hemodynamic delay (sec)"
  echo "  -sdf      file   : slice delay file"
  echo "  -dispersion tau  : hemodynamic dispersion (sec), -1 for sinusoid"
  echo "  -gammastem stem  : stem of volume with gamma function parameters"
  echo "  -snr       snr   : signal-to-noise ratio (-1 for all signal)"
  echo "  -pctsigch  pct   : percent signal change (wrt offset)"
  echo "  -scale     val   : overall scale"
  echo "  -nrows     n     : ($nrows)"
  echo "  -ncols     n     : ($ncols)"
  echo "  -ntp       n     : number time points"
  echo "  -roi r1 c1 r2 c2 : bounding box for activity"
  echo "  -firstslice n    : first anatomical slice"
  echo "  -nslices    n    : number of anatomical slices"
  echo "  -sfa   sfastem   : stem of selfreqavg output"
  echo "  -seed  seedval   : seed for random number generator (auto)"
  echo "  -monly mfile     : just create a matlab file"
  echo "  -umask umask     : set unix file permission mask"
  echo "  -version         : print version and exit"
exit 1;