#!/usr/bin/tcsh -f
# xcerebralseg

set VERSION = 'xcerebralseg mockbuild-local';

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 = 0;
set UseSamseg = 1
set ForceUpdate = 0
set SamsegRevert = 0
set ctab = $FREESURFER_HOME/FreeSurferColorLUT.txt 

# Seghead parameters
set thresh = 35
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 threads = $OMP_NUM_THREADS
set thresh1 = $thresh
set thresh2 = $thresh

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
  if($UseSamseg == 0) then
    set ud = `UpdateNeeded $seg1 $srcvolpath $xformpath`
    if($ud || $ForceUpdate) then
      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
  else
    set samsegoutdir = $SUBJECTS_DIR/$subject/mri/xcerebral.samseg
    set samsegseg = $samsegoutdir/seg.mgz
    set ud = `UpdateNeeded $samsegseg $srcvolpath`
    if($ud || $ForceUpdate) then
      set cmd = (fs_time samseg --i $srcvolpath --o $samsegoutdir --threads $threads --cpvcw)
      if($ForceUpdate) set cmd = ($cmd --force)
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) goto error_exit;
    endif
    set ud = `UpdateNeeded $seg1 $samsegseg`
    if($ud || $ForceUpdate) then
      # the samseg seg will have a bunch of segments that are not in the ca_label
      # output, so make sure to convert them. 
      set cmd = (mri_binarize --i $samsegseg --o $seg1 \
        --replaceonly 184 172  --replaceonly 183 172  --replaceonly 11300 3 \
        --replaceonly 12300 42 --replaceonly 915 165  --replaceonly 916 165 \
        --replaceonly 911 258  --replaceonly 907 258 --replaceonly 66 41 \
        --replaceonly 930 258  --replaceonly 259 258  --replaceonly 24 257 \
        --replaceonly 909 258  --replaceonly 126 16    --replaceonly 34 2 )
      set cmd = ($cmd --replaceonly 262 130); # Sinus becomes Air
      set cmd = ($cmd --replaceonly 908 258); # Eye muscles
      # replace CC with left WM so that there are no residual 192 voxels
      set cmd = ($cmd --replaceonly 192 2)
      # replace 267 PonsBellyArea with 174 Pons so as not to confuse gtmseg
      set cmd = ($cmd --replaceonly 267 174)
      if($SamsegRevert) then
        # Up to this point, the new seg will have the samseg set of labels as
        # the GCA seg. If SamsegRevert is set, then these last two will be replaced
        # as well making the set of labels the same in both
        set cmd = ($cmd --replaceonly 902 258); # Artery
        set cmd = ($cmd --replaceonly 914 257); # Vein
      endif
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) goto error_exit;
    endif
  endif
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 \
  --replaceonly 36 257 --replaceonly 66 257 --replaceonly 11300 257 --replaceonly 12300 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 267 --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 183 184 --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 --ctab $ctab --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 $ctab \
     --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 "--thresh":
      if($#argv < 1) goto arg1err;
      set thresh = $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 "--stats":
      set DoStats = 1;
      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 "--samseg":
      set UseSamseg = 1
      set srcvol = orig.mgz # nu can have holes in it
      breaksw
    case "--nosamseg":
    case "--no-samseg":
      set UseSamseg = 0
      # Note that if --samseg was used, this will not undo srcvol
      breaksw

    case "--samseg-revert":
      set SamsegRevert = 1
      breaksw

    case "--force":
     set ForceUpdate = 1
     breaksw
    case "--no-force":
     set ForceUpdate = 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($UseSamseg == 0 && $#seg1name == 0) then
  set xformpath = $SUBJECTS_DIR/$subject/mri/transforms/$xform
  if(! -e $xformpath) then
    echo "ERROR: cannot find $xformpath"
    echo " Try using the --samseg flag"
    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 : relative to subject/mri (default: $outvol)"
  echo " --atlas headgca (default: $atlas)"
  echo " --samseg : use samseg instead of GCA to provide the xcer seg"
  echo " --samseg-revert : make labels of samseg seg same as that of the GCA"
  echo " --m mergevol : merge with mergevol (default is $mergevol)"
  echo " --srcvol source intensity volume (default is $srcvol)"
  echo " --stats : run mri_segstats"
  echo " --seg1 seg1name : full-head seg (usually computed with mri_ca_label)"
  echo " --thresh thresh : threshold passed to mri_seghead ($thresh)"
  echo ""
  echo " --no-pons"
  echo " --no-vermis"
  echo " --threads nthreads : set OPENMP threads"
  echo " --force : force rerunning and overwriting results"

  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 

