#! /bin/csh -f

#
# bgmask-sess
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: nicks $
#    $Date: 2007/01/09 22:41:16 $
#    $Revision: 1.9 $
#
# 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: bgmask-sess,v 1.9 2007/01/09 22:41:16 nicks Exp $';

set cmdargs = ($argv);

set omask     = bgmask;
set fsd       = bold;
set funcstem  = f;
set epi       = 1;
set RunListFile = ();
set PerRun    = 0;
set PolyFitOrder = 2;
set noepi = 0;
set ExtReg    = ();
set nolog     = 0;
set PrintHelp = 0;
set monly     = 0;
set sumfile   = ();
set MLF       = ();
set infmt  = ();
set outfmt = nii.gz
if($?FSF_OUTPUT_FORMAT) then
  set outfmt = $FSF_OUTPUT_FORMAT;
endif

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

set SessList = `getsesspath $argv`;
if($status || $#SessList == 0) then
  getsesspath $argv 
  exit 1;
endif

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

if(! $nolog) then
  set LF = log/bgmask-sess.log
  if(-e $LF) mv $LF $LF.bak
else
  set LF = /dev/null
endif

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
uname -a  | tee -a $LF

if(! $monly) set MLF = /tmp/bgmask-sess.$$.m
rm -f $MLF

if($#sumfile) then 
  rm -f $sumfile
  echo "# $VERSION " >>$sumfile
  echo "# `date`" >>$sumfile
  echo "# PolyFitOrder $PolyFitOrder" >>$sumfile
  echo "# ExtReg $ExtReg" >>$sumfile
  echo "# FuncStem $funcstem" >>$sumfile
  echo "# " >>$sumfile
  echo "#    sessid     runid nbgmask      bgmn bgstd bgstdexp     fgmn fgstd fgstdw     pctnw" >> $sumfile
endif

#------------------------------------------------------------------#
set StartTime = `date`;
foreach sess ($SessList)
  set sessid = `basename $sess`;

  set RunList = `getrunlist $sess/$fsd $RunListFile`;
  if($status) then
    echo $RunList
    exit 1;
  endif

  if(! $PerRun) set RunList = $RunList[1];

  foreach Run ($RunList)

    set funcstem0 = $sess/$fsd/$Run/$funcstem
    set infmtuse = `getformat_from_stem $funcstem0 $infmt`
    if($status) then
      echo "$infmtuse" |& tee -a $LF
      exit 1;
    endif
    echo "Detected format as $infmtuse" |& tee -a $LF
    set func = $sess/$fsd/$Run/$funcstem.$infmtuse

    if(! $PerRun) then
      set maskdir = $sess/$fsd/masks
    else
      set maskdir = $sess/$fsd/$Run/masks
    endif
    mkdir -p $maskdir
    set outvol = $maskdir/$omask.$outfmt
  
    set extregstem = ();
    if($#ExtReg) then
      set extregstem = $sess/$fsd/$Run/$ExtReg
      if(! -e $extregstem"_000.hdr") then
        echo ERROR: cannot find $extregstem"_000.hdr";
        exit 1;
      endif
    endif

#--------------------------#
tee -a $MLF > /dev/null <<EOF
clear all;
tic;

QuitOnError = ~[$monly];
fspec = '$func';
extregstem = '$extregstem';
polyfitorder = [$PolyFitOrder];
noepi = $noepi;
bgmaskspec = '$outvol';
sessid = '$sessid';
runid  = '$Run';
sumfile = '$sumfile';

fprintf('%s Loading functional %s\n',sessid,runid);
fprintf('%s\n',fspec);
f = MRIread(fspec);
if(isempty(f))
  fprintf('ERROR: loading %s\n',fspec);
  if(QuitOnError) 
    quit; 
  end
  return;
end
f.volmat = fast_vol2mat(f.vol);

X = [];
Xpoly = fast_polytrendmtx(1,f.nframes,1,polyfitorder);
X = [X Xpoly];

if(~isempty(extregstem))
  extreg = fmri_ldbvolume(extregstem);
  if(isempty(extreg))
    fprintf('ERROR: could not load %s\n',extregstem);
    if(QuitOnError) quit; end
    return;
  end
  if(size(extreg,3)~=1) extreg = squeeze(extreg)'; %'
  else                  extreg = squeeze(extreg);
  end
  % Remove mean of External Regressor %
  extreg = extreg(:,1:3);
  extreg = extreg - repmat(mean(extreg), [f.nframes 1]);
  extreg = extreg./repmat(std(extreg), [f.nframes 1]);
  X = [X extreg];
end

if(0)
  % This has very little effect, so dont do it
  fprintf('SVD\n');
  tmpmn = mean(f.volmat);
  indtmp = find(tmpmn > 2*mean(tmpmn));
  fprintf('nSVD: %d\n',length(indtmp));
  ytmp = f.volmat(:,indtmp);
  ytmp = ytmp - X*(inv(X'*X)*X'*ytmp);
  [u s v] = fast_svd(ytmp);
  X = [X u(:,1:5)];
  ds = diag(s);
  ds = 100*ds/sum(ds);
  fprintf('SingVals: ');
  fprintf('%4.1f ',ds(1:5));
  fprintf('\n');
end

fprintf('  Detrending \n');
[beta, rvar] = fast_glmfitw(f.volmat,X);
fstd = sqrt(rvar);

fmn = fast_mat2vol(beta(1,:),f.volsize);
fstd = fast_mat2vol(fstd,f.volsize);

fprintf('  Computing background mask \n');
bgmask = f;
[bgmask.vol bgstd fstdfg fstdfgexp bgthresh fmnfg bgmB] = fast_bgmask(fmn,fstd,noepi);
indbg = find(bgmask.vol);
nbgmask = length(indbg);
bgmn = mean(fmn(indbg));
bgstdexp = sqrt((4/pi-1)*bgmn.^2); % Expected based on the mean

% Stocker's method
indbgB  = find(bgmB);
bgmnB   = mean(fmn(indbgB));
bgstdB  = mean(fstd(indbgB));
bgstdBexp = sqrt((4/pi-1)*bgmnB.^2); % Expected based on the mean

fprintf('%10s %s  n = %5d  bgmn = %8.3f bgstdexp = %8.3f  bgstd = %8.3f  fgstdexp = %8.3f  fgstd = %8.3f fgmn = %g\n',...
        sessid,runid,nbgmask,bgmn, bgstdexp, bgstd, fstdfgexp, fstdfg, fmnfg);

% percentage of variance due to non-white noise
pctnw = 100*(fstdfg.^2 - fstdfgexp.^2)/(fstdfg.^2);

localsumfile = sprintf('%s.bgsum',bgmaskspec);
fp = fopen(localsumfile,'w');
fprintf(fp,'sessid    %s\n',sessid);
fprintf(fp,'runid     %s\n',runid);
fprintf(fp,'fspec     %s\n',basename(fspec));
if(~isempty(extregstem))
  fprintf(fp,'extreg    %s\n',basename(extregstem));
end
fprintf(fp,'noepi     %d\n',noepi);
fprintf(fp,'poly      %d\n',polyfitorder);
fprintf(fp,'nbg       %d\n',nbgmask);
fprintf(fp,'bgmn      %f\n',bgmn);
fprintf(fp,'bgstd     %f\n',bgstd);
fprintf(fp,'bgstdexp  %f\n',bgstdexp);
fprintf(fp,'fgmn      %f\n',fmnfg);
fprintf(fp,'fgstd     %f\n',fstdfg);
fprintf(fp,'fgstdexp  %f\n',fstdfgexp);
fprintf(fp,'fgpctnw   %f\n',pctnw);
fprintf(fp,'bgmnB     %f\n',bgmnB);
fprintf(fp,'bgstdB    %f\n',bgstdB);
fprintf(fp,'bgstdBexp %f\n',bgstdBexp);
fclose(fp);

if(~isempty(sumfile))
  fp = fopen(sumfile,'a');
  fprintf(fp,'%15s %s  %5d   %8.3f %8.3f %8.3f   %8.3f %8.3f %8.3f    %4.1f\n',...
          sessid,runid,nbgmask,...
          bgmn,bgstd,bgstdexp, ...
          fmnfg,fstdfg,fstdfgexp, pctnw);
  fclose(fp);
end

if(0)
  fbg = f.volmat(:,indbg);
  x = [0:35];
  [h x] = hist(fbg(:),x);
  pdfest = h/trapz(x,h);
  fbgmn  = mean(fbg(:));
  fbgstd = std(fbg(:));
  rpdf = pdf_rayleigh(x,fbgmn);
  plot(x,pdfest,x,rpdf);
end

fprintf('  Saving \n');
MRIwrite(bgmask,bgmaskspec);

fprintf('  Done %g \n',toc);

EOF
#--------------------------#

  end # Runs

end # Session

set okfile = /tmp/bgmask-sess.$$.ok
rm -f $okfile
echo "fmri_touch('$okfile');" >> $MLF

if(! $monly) then
  cat $MLF | matlab -display iconic | tee -a $LF
  rm -f $MLF
  if(! -e $okfile) then
    echo "ERROR: in matlab execution"
    exit 1;
  endif
  rm -f $okfile
endif

echo "Started at $StartTime" | tee -a $LF
echo "Ended   at `date`"     | tee -a $LF
echo "bgmask-sess Done"    | tee -a $LF

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


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

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

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

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

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

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

    case "-mcextreg":
      set ExtReg = mcextreg;
      breaksw

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

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

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

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

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

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

    case "-perrun":
      set PerRun = 1;
      breaksw

    case "-epi":
      set noepi = 0;
      breaksw

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

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

    case "-g":
    case "-s":
    case "-sf":
    case "-d":
    case "-df":
      shift;
      # ignore getsesspath arguments 
      breaksw

    case "-cwd":
      # ignore getsesspath arguments 
      breaksw

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

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

end

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

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

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

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

############--------------##################
usage_exit:
  echo ""
  echo "USAGE: bgmask-sess"
  echo ""
  echo "Optional Arguments:";
  echo "   -o    outmask : use outmask instead of bgmask"
  echo "   -func funcstem : default is f"
  echo "   -fsd  fsd : use fsd instead of bold"
  echo "   -infmt fmt : input format <auto>, nii, mgh, bhdr, etc"
  echo "   -outfmt fmt : output format <nii.gz>, mgh, bhdr, etc"
  echo "   -rlf  rlf : run list file"
  echo "   -extreg stem : external regressors"
  echo "   -mcextreg  : same as -extreg mcextreg"
  echo "   -perrun   : analyze each run separately (RRR/masks)"
  echo "   -noepi    : do not assume EPI"
  echo "   -sum sumfile : save summary of all runs in sumfile"
  echo ""
  echo "Session Arguments (Required)"
  echo "   -sf sessidfile  ..."
  echo "   -df srchdirfile ..."
  echo "   -s  sessid      ..."
  echo "   -d  srchdir     ..."
  echo ""
  echo "Session Arguments (Optional)"
  echo "   -version       : print version and exit"
  echo "   -debug"
  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