#! /bin/csh -f

#
# mergecontrasts-sess - AND conjunction of two contrasts
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: nicks $
#    $Date: 2007/01/09 22:41:18 $
#    $Revision: 1.6 $
#
# 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: mergecontrasts-xsess,v 1.6 2007/01/09 22:41:18 nicks Exp $';

set inputargs = ($argv);

set mergedcontrast = ();
set conjunction = ();
set map = sig;
set analysis    = ();
set space = ();
set spacedir = ();
set hemi = ();
set hemisuf = ();
set isxavg       = ();
set sesslist = ();
set outsess = ();

set ncontrasts   = 0;
set contrastlist = ();
set sessidlist   = ();
set threshlist   = ();
set taillist     = ();
set complist     = ();

set monly = 0;
set MLF = ();
set nolog = 0;
set checkvolumes = 0;

set PrintHelp = 0;

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;
endif

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

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

## Get functional subdirectory from the info file ##
set infofile = $analysis/analysis.info
if(! -e $infofile) then
  echo "ERROR: cannot find $infofile"
  exit 1;
endif
set fsd = `cat $analysis/analysis.info | awk '{if($1 == "fsd") print $2}'`;

## Look for nolog option ##
set n = `echo $inputargs | grep nolog | wc -l` 
if($n != 0) set nolog = 1;

##### Create a log file ######
if(! $nolog) then
  set logdir = `pwd`/log;
  mkdir -p $logdir
  if(! -e $logdir) then
    echo "ERROR: could not create $logdir"
    exit 1;
  endif
  set LF = $logdir/mergecontrasts-xsess.log
  if(-e $LF) mv $LF $LF.old
else
  echo "No log file"
  set LF = /dev/null
endif

echo "mergecontrasts-xsess Logfile is $LF"

echo "mergecontrasts-xsess log file" >> $LF
echo $VERSION >> $LF
id            >> $LF
pwd           >> $LF
echo "$0"     >> $LF
echo $inputargs  >> $LF
uname -a      >> $LF
date          >> $LF

set inbhdr  = ();
set instems = ();
@ nthcontrast = 1;
foreach contrast ($contrastlist)
  set sess = $sesslist[$nthcontrast];
  set indir = $sess/$fsd/$analysis/$spacedir-$isxavg/$contrast
  if(! -d $indir) then
    echo "ERROR: $indir does not exist"   |& tee -a $LF
    exit 1;
  endif
  set instemtmp = $indir/$map$hemisuf;
  set instems = ($instems $instemtmp);
  if(-e $instemtmp.bhdr) set inbhdr = $instemtmp.bhdr
  echo $instemtmp
  @ nthcontrast = $nthcontrast + 1;
end

if($checkvolumes) then
  foreach invol ($instems)
    echo "Checking $invol"
    bfvcheck $invol
    if($status) then
      echo "ERROR: there appears to be a problem with input volume $instem1"
      exit 1;
    endif
  end
endif

set outdir = $outsess/$fsd/$analysis/$spacedir-$isxavg/$mergedcontrast
mkdir -p $outdir
if(! -d $outdir) then
  echo "ERROR: cannot create $outdir"|& tee -a $LF
  exit 1;
endif
set outstem = $outdir/$map$hemisuf
set ovlpfile = $outstem.ovlp

if(! $monly ) set MLF = $outdir/mergcon-$$.m
rm -f $MLF  

echo $instems

#----------------------------------------------------------#
tee $MLF > /dev/null <<EOF

  monly = $monly;
  if(monly) QuitOnError = 0;
  else      QuitOnError = 1;
  end

  outstem = '$outstem';
  conjunction = '$conjunction';
  ovlpfile = '$ovlpfile';

  instems = splitstring('$instems');
  threshs = [$threshlist];
  tails   = splitstring('$taillist');
  comps   = [$complist];

  fprintf('\n');
  nover = [];
  ninputs = size(instems,1);
  for n = 1:ninputs
    instem = deblank(instems(n,:));
    fprintf('Loading %s\n',instem);
    m = fmri_ldbvolume(instem);

    % check the size
    if(n==1) 
      sz0 = size(m);
    else
      sz = size(m);
      for i = 1:3
        if(sz(i) ~= sz0(i))
         fprintf('ERROR: volume dimensions do not match\n');
	 fprintf('  %s\n',deblank(instems(1,:)));
	 fprintf('  %s\n',deblank(instems(n,:)));
        end
      end
    end

    if(comps(n)) m = sign(m).*(1-10.^(-abs(m))); end

    % mthresh is a map of the nth input whose voxels meet thresh %
    mthresh = fast_voxthresh(m,threshs(n),deblank(tails(n,:)));
    novertmp = length(find(mthresh > 0.5));
    nover = [nover novertmp ];
    fprintf('  Found %d voxels\n',novertmp);

    if(n==1) 
      mmerged = mthresh;
    else
      switch(conjunction)
      case 'and',   mmerged = mthresh & mmerged;
      case 'or',    mmerged = mthresh | mmerged;
      case 'andor', 
        % This only makes sense for 2 inputs %
        map_both = mthresh & mmerged;
	ind_both = find(map_both);
        indA = find(mmerged);
        indB = find(mthresh);
        mmerged_tmp = zeros(size(mmerged));
        mmerged_tmp(indA)     = +1;
        mmerged_tmp(indB)     = -1;
        mmerged_tmp(ind_both) = +2;
        mmerged = mmerged_tmp;
      end
    end

  end
  fprintf('Saving result to %s\n',outstem);
  fmri_svbvolume(mmerged,outstem);

  noverlap = length(find(mmerged > 0.5));
  fprintf('nfinal = %d\n',noverlap);

  fid = fopen(ovlpfile,'w');
  if(fid == -1) 
    fprintf('ERROR: could not open %s\n',ovlpfile);
  else
    fprintf(fid,'%d ',[nover noverlap]);
    fprintf(fid,'\n');
    fclose(fid);
  end

  if(~monly) 
    fprintf('quitting matlab\n');
    quit; 
  end

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

cat $MLF >> $LF

if(! $monly ) then
    echo "Starting matlab"
    cat $MLF | matlab -display inconic |& tee -a $LF
    rm -f $MLF  
    if($#inbhdr != 0) cp $inbhdr $outstem.bhdr
endif



echo mergecontrasts-xsess Done |& tee -a $LF

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

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

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

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

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

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

    case "-contrast":
    case "-c":
      if( $#argv < 5) then
        echo "ERROR: -contrast requires 5 arguments:"
        echo "        contrastname, sessid, threshold, tail, and complement"
        exit;
      endif
      set contrastlist = ($contrastlist $argv[1]); shift;
      set sesslist     = ($sesslist     $argv[1]); shift;
      set threshlist   = ($threshlist   $argv[1]); shift;
      set taillist     = ($taillist     $argv[1]); shift;
      set complist     = ($complist     $argv[1]); shift;
      breaksw

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

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

    case "-space":
      if ( $#argv == 0) goto arg1err;
      set space = $argv[1]; shift;
      if($space != tal && $space != sph) then
        echo "ERROR: space must be either tal or sph"
        exit 1;
      endif
      breaksw

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

    case "-hemi":
      if ( $#argv == 0) goto arg1err;
      set hemi = $argv[1]; shift;
      if($hemi != lh && $hemi != rh) then
        echo "ERROR: hemi must be either lh or rh"
        exit 1;
      endif
      breaksw

    case "-isxavg"
      if ( $#argv == 0) goto arg1err;
      set isxavg = $argv[1]; shift;
      if("$isxavg" != random & "$isxavg" != fixed) then
        echo "ERROR: isxavg must be either random or fixed"
        exit 1;
      endif
      if("$isxavg" == random) set isxavg = rfx;
      if("$isxavg" == fixed)  set isxavg = ffx;
      breaksw

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

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

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

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

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

    case "-nolog"
    case "-cwd"
      breaksw

    case "-s":
    case "-sf":
    case "-df":
    case "-d":
    case "-g":
      shift;
      breaksw

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

end

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

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

############--------------##################
arg4err:
  echo "ERROR: flag $flag requires four arguments"
  exit 1
############--------------##################

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

  if($#analysis == 0) then
    echo "ERROR: no analysis specified"
    exit 1
  endif

  if($#map == 0) then
    echo "ERROR: no map specified"
    exit 1
  endif

  set ncontrasts = $#contrastlist 
  if($ncontrasts < 2) then
    echo "ERROR: must specify at least 2 contrasts"
    exit 1
  endif

  if($#mergedcontrast == 0) then
    echo "ERROR: must specify merged contrast name "
    exit 1
  endif

  if($#conjunction == 0) then
    echo "ERROR: must specify conjunction (and, or, andor) "
    exit 1
  endif

  if("$conjunction" != "and" && "$conjunction" != "or" && "$conjunction" != "andor") then
    echo "ERROR: conjunction = $conjunction, must be 'and', 'or', or 'andor'"
    exit 1;
  endif

  foreach tail ($taillist)
    if($tail != "pos" && $tail != "neg" && $tail != "abs") then
      echo "ERROR: tail = $tail, must be pos, neg, or abs"
      exit 1;
    endif
  end

  foreach comp ($complist)
    if($comp != 1 && $comp != 0) then
      echo "ERROR: complement = $comp, must be 0 or 1"
      exit 1;
    endif
  end

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

  if($#hemi == 0 && $space == "sph") then
    echo "ERROR: for space sph, need -hemi "
    exit 1;
  endif

  if($#hemi != 0 && $space != "sph") then
    echo "ERROR: for -hemi, space must be sph"
    exit 1;
  endif

  if($#hemi != 0) set hemisuf = "-$hemi";
  if($#space != 0 && $#spacedir == 0) set spacedir = $space;

  if($#outsess == 0)  set outsess = $sesslist[1];
  echo "INFO: output session is $outsess"

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

############--------------##################
usage_exit:
  echo "USAGE: mergecontrasts-xsess"
  echo "Options:";
  echo "   -mergedcontrast name of new contrast"
  echo "   -conjunction name : and or or "
  echo "   -analysis analysisname : name of session-level functional analysis";
  echo "   -contrast   name sessid thresh tail complement"
  echo "   -map       mapname"
  echo "   -space  spacename : tal or sph"
  echo "   -spacedir spacedir : if not tal or sph"
  echo "   -isxavg ffx or rfx"
  echo "   -outsess output sess (default is first)"
  echo "   -version       : print version and exit"
  echo "   -help          : print help and exit"
  echo "   -umask umask   : set unix file permission mask"

  if($PrintHelp) \
  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

This progam will merge two or more contrasts across sessions. The
merging is done by comparing each thresholded contrast on a voxel-by-
voxel basis. The result depends upon the conjunction chosen. For the
'and' conjunction, each contrast must meet the threshold criteria at a
voxel for the voxel to be set to 1.  For the 'or' conjunction, the
voxel is set to 1 if any contrast meets the threshold
criteria. Otherwise, the voxel is set to zero. 

The output is a map with the same name as the input map which will be
placed in a new contrast directory in the given output session. If no
output session is specified then the first session in the contrast
list is used. The merged contrast can be accessed in the same was as
any other contrast.  However, remember that the voxel values of the
merged contrast will be either 0 or 1, so set thresholds
appropriately.

The contrast and threshold criteria are set with the -contrast flag.
The first argument is the name of the contrast. The second is the name
of the session, The third is the threshold to apply. If the map is a
significance map, then set the threshold based on log10(p) (eg, 2 =
.01). The fourth argument is the tail, which can be either pos, neg,
or abs for positive, negative, or ignore sign. The final arugment
allows the user to use a complement. Ie, if the fourth argument is
non-zero and the threshold criteria is met, then a 0 is input to the
conjunction.

There is a .ovlp file created with the following three numbers:
 (a) number of voxels in contrast 1 that met threshold criteria
 (b) number of voxels in contrast 2 that met threshold criteria
 (c) number of voxels that overlapped

Example:

mergecontrasts-xsess \
  -mergedcontrast conmerge-or \
  -a sirp-mcsm5x \
  -conj or \
  -c t2va grp-norm  2 abs 0 \
  -c t5va grp-schiz 2 abs 0 \
  -map sig -space tal -isxavg fixed

This will create a new contrast called conmerge-or by oring the
t2va contrast from session/group grp-norm and t5va from grp-schiz.
The result will be stored in grp-norm (but this could have been
changed with -outsess.

Notes:

1. There is no -s or -sf input because the sessions are specified
as arguments to the -c option.

2. Currently does not accept -d or -df options. This means that
the sessions must be subdirectories of the study directory. I can
change this once I have time.