#!/bin/tcsh -f
# xcerebralseg

set VERSION = 'xcerebralseg dev';

set outvol = apas+head.mgz
set subject = ();
#set atlas = $FREESURFER_HOME/average/aseg+spmhead.ixi.gca
set atlas = $FREESURFER_HOME/average/aseg+spmhead+vermis+pons.ixi.gca
set srcvol = nu.mgz
set xform = talairach.m3z
set mergevol = aparc+aseg.mgz
set DoPons = 1;
set DoVermis = 1;
set seg1name = ()
set DoStats = 1;

# Seghead parameters
set thresh1 = 20;
set thresh2 = 20;
set nhitsmin = 2;

if ($?FS_OMP_NUM_THREADS) then
  setenv OMP_NUM_THREADS $FS_OMP_NUM_THREADS
else
  setenv OMP_NUM_THREADS 1
endif

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

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

set StartTime = `date`;
set tSecStart = `date '+%s'`;
set year  = `date +%Y`
set month = `date +%m`
set day   = `date +%d`
set hour   = `date +%H`
set min    = `date +%M`

set sd = $SUBJECTS_DIR/$subject

if($#tmpdir == 0) then
  set tmpdir = `fs_temp_dir --scratch`
endif
mkdir -p $tmpdir

# Set up log file
if($#LF == 0) set LF = $sd/scripts/xcerebralseg.log
if($LF != /dev/null) rm -f $LF
echo "Log file for xcerebralseg" >> $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
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
uname -a  | tee -a $LF

#========================================================
set outseg = $SUBJECTS_DIR/$subject/mri/$outvol
set segbase = `fname2stem $outvol`

# Perform the labeling
if($#seg1name == 0) then
  set seg1 = $tmpdir/seg1.mgh
  set cmd = (mri_ca_label -align -nobigventricles \
     $srcvolpath $xformpath $atlas $seg1)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
endif

# Replace all cortex GM and cerebral WM with xCSF. It wont matter
# that stuff in the brain is replaced because that will get replaced
# in the merge anyway
set seg2 = $tmpdir/seg2.mgh
set cmd = (mri_binarize --i $seg1 --o $seg2 \
  --replaceonly 3 257 --replaceonly 42 257 --replaceonly 2 257 --replaceonly 41 257)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) goto error_exit;

# Merge the FS brain seg with this seg
set seg3 = $tmpdir/seg3.mgh
set cmd = (mergeseg --src $seg2 --merge $mergevolpath --o $seg3)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) goto error_exit;

if($DoPons) then
  # Make a mask of pons
  set pons = $tmpdir/pons.mgh
  set cmd = (mri_binarize --i $seg1 --match 174 --o $pons --binval 174)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  # Merge pons
  set seg3b = $tmpdir/seg3b.mgh
  set cmd = (mergeseg --src $seg3 --merge $pons --o $seg3b)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  set seg3 = $seg3b;
endif

if($DoVermis) then
  # Make a mask of vermis
  set vermis = $tmpdir/vermis.mgh
  set cmd = (mri_binarize --i $seg1 --match 172 --o $vermis --binval 172)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  # Merge pons
  set seg3c = $tmpdir/seg3c.mgh
  set cmd = (mergeseg --src $seg3 --merge $vermis --o $seg3c)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  set seg3 = $seg3c;
endif

# Create a seg of the head for filling in voxels and masking stuff out
set seghead = $tmpdir/seghead.mgz
set cmd = (mri_seghead --invol $srcvolpath --outvol $seghead \
  --fill 1 --thresh1  $thresh1 --thresh2  $thresh2 --nhitsmin $nhitsmin);
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) exit 1;

# Replace unlabeled voxels in the head with head seg
set seg4 = $tmpdir/seg4.mgh
set cmd = (mri_binarize --i $seg3 --o $seg4 --replaceonly 0 258 --mask $seghead);
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) goto error_exit;

# Replace unlabeled voxels in the head with head seg
set cmd = (mri_mask $seg4 $seghead $outseg)
echo $cmd | tee -a $LF
$cmd | tee -a $LF
if($status) goto error_exit;

# Run segstats to get stats of the extra cerebral structures.  This is
# done for informational purposes only.  Values for other structures
# are there for reference.
if($DoStats) then
  set sdir = $SUBJECTS_DIR/$subject
  set sum =  $sdir/stats/$segbase.stats 
  set cmd = (mri_segstats --seg $outseg --sum $sum \
     --pv $sdir/mri/nu.mgz --empty --brainmask $sdir/mri/brainmask.mgz \
     --in-intensity-name nu \
     --in $sdir/mri/nu.mgz  --in-intensity-units MR --etiv  \
     --ctab $FREESURFER_HOME/FreeSurferColorLUT.txt \
     --subject $subject --id 2 41 4 5 14 15 43 44 130 165 257 258)
  if($DoVermis) set cmd = ($cmd --id 172)
  if($DoPons)   set cmd = ($cmd --id 174)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
endif

#========================================================

# Cleanup
if($cleanup) rm -rf $tmpdir

# Done
echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "Xcerebralseg-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Xcerebralseg-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "xcerebralseg Done" |& tee -a $LF
exit 0

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

############--------------##################
error_exit:
echo "ERROR:"

exit 1;
###############################################

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

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

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

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

    case "--threads":
    case "--nthreads":
      setenv OMP_NUM_THREADS $argv[1]; shift;
      breaksw

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

    case "--no-stats":
      set DoStats = 0;
      breaksw

    case "--nopons":
    case "--no-pons":
      set DoPons = 0;
      breaksw

    case "--novermis":
    case "--no-vermis":
      set DoVermis = 0;
      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 "--tmp":
    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($#subject == 0) then
  echo "ERROR: must spec subject"
  exit 1;
endif
if(! -e $SUBJECTS_DIR/$subject) then
  echo "ERROR: cannot find $subject"
  exit 1;
endif

if($#seg1name) then
  set seg1 = $SUBJECTS_DIR/$subject/mri/$seg1name
  if(! -e $seg1) then
    echo "ERROR: cannot find $seg1"
    exit 1;
  endif
endif

if($#seg1name == 0) then
  set xformpath = $SUBJECTS_DIR/$subject/mri/transforms/$xform
  if(! -e $xformpath) then
    echo "ERROR: cannot find $xformpath"
    exit 1;
  endif
endif

set srcvolpath = $SUBJECTS_DIR/$subject/mri/$srcvol
if(! -e $srcvolpath) then
  echo "ERROR: cannot find $srcvolpath"
  exit 1;
endif

set mergevolpath = $SUBJECTS_DIR/$subject/mri/$mergevol
if(! -e $mergevolpath) then
  echo "ERROR: cannot find $mergevolpath"
  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 "xcerebralseg --help"
  echo " --s subject (required)"
  echo " --o outvol (default: $outvol)"
  echo " --atlas headgca (default: $atlas)"
  echo " --m mergevol : merge with mergevol (default is $mergevol)"
  echo " --srcvol source intensity volume (default is $srcvol)"
  echo " --no-stats : do not  run mri_segstats"
  echo " --seg1 seg1name : full-head seg (usually computed with mri_ca_label)"
  echo ""
  echo " --no-pons"
  echo " --no-vermis"
  echo " --threads nthreads : set OPENMP threads"
  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

Performs labeling of extracerebral structures including sulcal CSF,
skull/bone, head soft tissue, and air inside the head (eg, sinuses).
This is merged with the aparc+aseg.mgz segmentation to give a whole
head segmentation.  It uses a GCA (aseg+spmhead.ixi.gca) generated
from 79 subjects from the IXI database (www.brain-development.org);
see below for a demographic break down of the 79. The 79 were analyzed
using SPM New Segment, and then FreeSurfer code was used to generate
the GCA based on the SPM segmentation. It has been tested against a CT
data set and performs about as well as SPM in terms of detecting the
skull. Still, the results are far from perfect, so one should still
treat them as approximate. This segmentation is primarily intended to
be used to create nuisance variables for fMRI and PET. Note that the
"Skull" segmentation includes any kind of bone.

ntotal 79 
males 34 
females 45 
white 39 
asian 14 
black 13 
chinese 13 
Philips-1.5T 39  (Guys)
Philips-3T 22    (HH)
GE-1.5T 18       (IOP)
mean age: 58.0995
min age: 20.21
max age: 86.32
age decade break down
 20-30  5 
 30-40 10 
 40-50 10 
 50-60 13 
 60-70 17 
 70-80 16 
 80-90  8