#!/bin/tcsh -f
# vlrmerge

set VERSION = '$Id: vlrmerge,v 1.2.2.2 2011/03/31 21:32:21 greve Exp $';

set outvol = ();
set invol = ();
set lhsurf = ();
set rhsurf = ();
set lrvol = ();
set subcortmask = ();
set DoCorrection = 0; # Correction for multiple comparsisons

set tmpdir = ();
set cleanup = 1;
set LF = ();

set inputargs = ($argv);
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:

set outdir = `dirname $outvol`

mkdir -p $outdir
pushd $outdir > /dev/null
set outdir = `pwd`;
popd > /dev/null

if($#tmpdir == 0) set tmpdir = $outdir/tmpdir.vlrmerge
mkdir -p $tmpdir

if($#LF == 0) then
  set stem = `fname2stem $outvol`
  set stem = `basename $stem`
  set LF = $outdir/vlrmerge.$stem.log
endif
if($LF != /dev/null) rm -f $LF

echo "Log file for vlrmerge" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
echo "" | tee -a $LF
echo $VERSION | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
uname -a  | tee -a $LF

# 1. Map the lh and rh into the volume
set surfvol = $tmpdir/junk.mgh
set cmd = (mri_surf2vol --fstal 2 \
  --template $SUBJECTS_DIR/fsaverage/mri.2mm/orig.mgz \
  --hemi lh --surfval $lhsurf --o $surfvol --fillribbon)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
set cmd = (mri_surf2vol --fstal 2 --merge $surfvol \
   --hemi rh --surfval $rhsurf  --o $surfvol --fillribbon)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
# Save copy of LRVol if desired
if($#lrvol) then
  set cmd = (mri_convert $surfvol $lrvol)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif
# 2. Invert the subcort mask
set subcortmaskinv = $tmpdir/subcortmaskinv.nii
set cmd = (mri_binarize --min 0.5 --i $subcortmask --inv \
  --o $subcortmaskinv)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
# 3. Mask the surf-in-vol to remove subcortex (will likely remove some cortex too)
set cmd = (mri_mask $surfvol $subcortmaskinv $surfvol)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
# 4. Mask the non-subcort out of the volume analysis
set involmasked = $tmpdir/involmasked.nii
set cmd = (mri_mask $invol $subcortmask $involmasked)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;
# 5. Add the subcort-only vol and the non-sub-cort surf-in-vol
set cmd = (mri_concat --sum $surfvol $involmasked --o $outvol)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;

if($DoCorrection) then
  # Bonferroni correct for multiple comparisons across the 3
  # spaces. Assumes that correction within a space have already been
  # done
  set cmd = (fscalc $outvol bcor 3 -o $outvol)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif

if($cleanup) rm -r $tmpdir

date | tee -a $LF
echo "vlrmerge done" | tee -a $LF

exit 0

###############################################

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

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

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

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

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

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

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

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

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

    case "--correct":
      # Bonferroni correction for multiple comparsisons across 3 spaces
      # Assumes that inputs have already been corrected
      set DoCorrection = 1; 
      breaksw

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

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      breaksw

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

    case "--nocleanup":
      set cleanup = 0;
      breaksw

    case "--cleanup":
      set cleanup = 1;
      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:

if($#invol == 0) then
  echo "ERROR: must specify input volume"
  exit 1;
endif
if($#lhsurf == 0) then
  echo "ERROR: must specify lh surf"
  exit 1;
endif
if($#rhsurf == 0) then
  echo "ERROR: must specify rh surf"
  exit 1;
endif
if($#outvol == 0) then
  echo "ERROR: must specify out volume"
  exit 1;
endif
if($#subcortmask == 0) then
  echo "ERROR: must specify a subcortical mask "
  exit 1;
endif


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

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

############--------------##################
usage_exit:
  echo ""
  echo "vlrmerge"
  echo "  --o outvol"
  echo "  --v vol"
  echo "  --lh lhsurf"
  echo "  --rh rhsurf"
  echo "  --subcortmask mask (--scm)"
  echo "  --correct : bonferroni correct across 3 spaces"
  echo "  "
  echo "  --lrvol lrvol : output a volume with only lh and rh"
  echo ""
  echo "  --help"
  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

This program merges volume and surface data into a volume for THE
PURPOSES OF DISPLAY ONLY. Do not attempt to draw inferences from this
image!  This program is meant to be used as part of the FSFAST group
fMRI analysis. In FSFAST, one analyzes data in three spaces: mni305,
left hemisphere, and right hemisphere, and group analysis results are
produced in each space. This program allows these maps to be merged
into a single volume using the surface-based analysis for the cortex
and volume-based analysis for subcortical. This can allow you to look
at all results inside the volume.

The usage would be something like:
1. Create and run FSFAST analysis in each space
2. Run isxconcat-sess for each space
3. Run mri_glmfit to create sig.nii in each space 

vlrmerge --o vlr.nii.gz \
  --v  analysis.mni305/contrast/glmdir/groupcontrast/sig.nii.gz \
  --lh analysis.lh/contrast/glmdir/groupcontrast/sig.nii.gz \
  --rh analysis.rh/contrast/glmdir/groupcontrast/sig.nii.gz \
  --subcortmask  analysis.mni305/subcort.mask.nii.gz

View with
  tkmedit fsaverage orig.mgz -aparc+aseg -ov vlr.mgh

If you do corrections for multiple comparisons for each space, you
can combine the cluster-corrected images into the volume and 
correct for comparisons across the three spaces with something like

vlrmerge --o vlr.corrected.mgh --correct \
  --vol analysis.mni305/contrast/glmdir/osgm/grf.th3.pos.sig.cluster.nii.gz \
  --lh  analysis.lh/contrast/glmdir/osgm/cache.th30.pos.sig.cluster.nii.gz  \
  --rh  analysis.rh/contrast/glmdir/osgm/cache.th30.pos.sig.cluster.nii.gz  \
  --scm analysis.mni305/subcort.mask.nii.gz)

