#! /bin/tcsh -f

#
# mkheadsurf
#
#
# Original Author: Doug Greve
#
# Copyright © 2021 The General Hospital Corporation (Boston, MA) "MGH"
#
# Terms and conditions for use, reproduction, distribution and contribution
# are found in the 'FreeSurfer Software License Agreement' contained
# in the file 'LICENSE' found in the FreeSurfer distribution, and here:
#
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
#
# Reporting: freesurfer@nmr.mgh.harvard.edu
#
#

# mkheadsurf

set VERSION = 'mkheadsurf dev';

set inputargs = ($argv);

umask 2;
set noseghead = 0;
set subjid = ();
set srcvol     = T1.mgz;
set headvol    = seghead.mgz;
set headsurf   = seghead;
set smheadsurf = smseghead;
set hemi = lh;
set fhi = ();
set thresh1 = 20;
set thresh2 = 20;
set nhitsmin = 2;
set fillval  = 1;
set NSmooth = 10;
set PrintHelp = 0;
set ndilate = 0;
set nerode = 0
set InputVol = ()
set OutputVol = ()
set OutputSurf = ()
set OrMask = ()
set LF = ()

set DoInflate = 0;
set DoCurv = 0
set Rescale = 1
set FillHolesIslands = 1
set UseMRITess = 0
set UseMarchingCubes = 1

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

set PWD = `getpwdcmd`;
if($status) exit 1;

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

mkdir -p $outvoldir $outsurfdir
if($#subjid) then
  set subjdir = $SUBJECTS_DIR/$subjid
  mkdir -p $subjdir/scripts
  if($#LF == 0) set LF = $subjdir/scripts/mkheadsurf.log
else
  if($#LF == 0) set LF = $outvoldir/mkheadsurf.log
endif
if(-e $LF) mv $LF $LF.old
echo "INFO: log file is $LF"

echo "Log file for mkheadsurf" >> $LF
date >> $LF
pwd >> $LF
echo $0 >> $LF
echo $inputargs >> $LF

set StartTime = `date`;

if(! $noseghead ) then
  #----- Segment the head ---------#
  set cmd = (mri_seghead );
  set cmd = ($cmd --invol    $InputVol)
  set cmd = ($cmd --outvol   $OutputVol)
  set cmd = ($cmd --fill     $fillval);
  set cmd = ($cmd --thresh1  $thresh1);
  set cmd = ($cmd --thresh2  $thresh2);
  set cmd = ($cmd --nhitsmin $nhitsmin);
  if($#fhi) set cmd = ($cmd --fhi $fhi);
  if($Rescale) set cmd = ($cmd --rescale)
  if(! $Rescale) set cmd = ($cmd --no-rescale)
  if($FillHolesIslands) set cmd = ($cmd --fill-holes-islands)
  if(! $FillHolesIslands) set cmd = ($cmd --no-fill-holes-islands)
  if($#OrMask) set cmd = ($cmd --or-mask $OrMask)
  echo "--------------------------------" |& tee -a $LF
  date |& tee -a $LF
  pwd |& tee -a $LF
  echo $cmd |& tee -a $LF
  echo "--------------------------------" |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) then
    echo "ERROR: mri_seghead" |& tee -a $LF
    exit 1;
  endif
  if($ndilate != 0 || $nerode != 0) then
    set cmd = (mri_binarize --i $OutputVol --min 0.5 \
      --o $OutputVol --binval $fillval)
    if($dilate > 0) set cmd = ($cmd --dilate $ndilate)
    if($nerode > 0) set cmd = ($cmd --erode $nerode)
    echo $cmd |& tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
else
  if(! -e $OutputVol) then
    echo "ERROR: $subjdir/mri/$headvol does not exist, and you "
    echo "have not chosen to create it"
    exit 1;
  endif
endif

#----- Tessellate Head ---------#
# Get full path because surface writing may try to put lh/rh in front of the file name
set OutputSurf = `getfullpath $OutputSurf`
if($UseMRITess) then
  set cmd = (mri_tessellate $OutputVol $fillval $OutputSurf)
endif
if($UseMarchingCubes) then
  set cmd = (mri_mc $OutputVol $fillval $OutputSurf)
endif

echo "--------------------------------" |& tee -a $LF
date |& tee -a $LF
pwd |& tee -a $LF
echo $cmd |& tee -a $LF
echo "--------------------------------" |& tee -a $LF
$cmd |& tee -a $LF
if($status) then
  echo "ERROR: mri_tessellate" |& tee -a $LF
  exit 1;
endif

#----- Smooth the tessellation ---------#
set cmd = (mris_smooth -n $NSmooth)
if($DoCurv) then
  set cmd = ($cmd -b $outsurfdir/area.$headsurf -c $outsurfdir/curv.$headsurf)
else
  set cmd = ($cmd -nw)
endif
set cmd = ($cmd $OutputSurf $OutputSurf)
echo "--------------------------------" |& tee -a $LF
date |& tee -a $LF
pwd |& tee -a $LF
echo $cmd |& tee -a $LF
echo "--------------------------------" |& tee -a $LF
$cmd |& tee -a $LF
if($status) then
  echo "ERROR: mris_smooth" |& tee -a $LF
  exit 1;
endif

#----- Inflate ---------#
if($DoInflate) then
  set cmd = (mris_inflate -n 10 -sulc $outsurfdir/sulc.$headsurf \
     $OutputSurf $OutputSurf.inflated)
  echo "--------------------------------" |& tee -a $LF
  date |& tee -a $LF
  pwd |& tee -a $LF
  echo $cmd |& tee -a $LF
  echo "--------------------------------" |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
endif

echo "Started at: $StartTime" |& tee -a $LF
echo "Ended   at: `date`"     |& tee -a $LF

echo "mkheadsurf done"  |& tee -a $LF


exit 0

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

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

  set flag = $argv[1]; shift;

  switch($flag)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    case "-rescale":
      set Rescale = 1
      breaksw
    case "-no-rescale":
      set Rescale = 0
      breaksw

    case "-no-fill-holes-islands":
      set FillHolesIslands = 0
      breaksw
    case "-fill-holes-islands":
      set FillHolesIslands = 1
      breaksw

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

    case "-curv":
      set DoCurv = 1
      breaksw
    case "-nocurv":
    case "-no-curv":
      set DoCurv = 0
      breaksw

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

    case "-tess":
      set UseMRITess = 1
      set UseMarchingCubes = 0
      breaksw
    case "-mc":
      set UseMRITess = 0
      set UseMarchingCubes = 1
      breaksw

    case "-inflate":
      set DoInflate = 1;
      breaksw

    case "-noinflate":
    case "-no-inflate":
      set DoInflate = 0;
      breaksw

    case "-log":
    case "--log":
      if($#argv < 1) goto arg1err;
      set LF = $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 "-g":
    case "-s":
    case "-sf":
    case "-d":
    case "-df":
      shift;
      # ignore getsesspath arguments
      breaksw

    case "-cwd":
      # ignore getsesspath arguments
      breaksw

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

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

end

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

############--------------##################
check_params:
  if($#subjid != 0) then
    if(! -e $SUBJECTS_DIR/$subjid ) then
      echo "ERROR: cannot find subject $subjid in $SUBJECTS_DIR"
      exit 1;
    endif
    set InputVol   = $SUBJECTS_DIR/$subjid/mri/$srcvol
    set OutputVol  = $SUBJECTS_DIR/$subjid/mri/$headvol
    set OutputSurf = $SUBJECTS_DIR/$subjid/surf/$hemi.$headsurf
  endif

  if($#InputVol == 0) then
    echo "ERROR: must spec input volume"
    exit 1;
  endif
  if($#OutputVol == 0) then
    echo "ERROR: must spec output volume"
    exit 1;
  endif
  if($#OutputSurf == 0) then
    echo "ERROR: must spec output surf"
    exit 1;
  endif
  if(! -e $InputVol ) then
    echo "ERROR: cannot find input volume $InputVol"
    exit 1;
  endif

  set outvoldir = `dirname $OutputVol`
  set outsurfdir = `dirname $OutputSurf`


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

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

############--------------##################
usage_exit:
  echo ""
  echo "mkheadsurf"
  echo ""
  echo "Required Arguments:";
  echo "  -i inputvol"
  echo "  -o outputvol"
  echo "  -surf ouputsurf"
  echo ""
  echo "  -s subjectid"
  echo "    Sets inputvol to be subject/mri/T1.mgz"
  echo "    Sets outputvol to be subject/mri/seghead.mgz"
  echo "    Sets outputsurf to be subject/surf/lh.smseghead"
  echo ""
  echo "Other Arguments (Optional)"
  echo "   -nsmooth N : number of smoothing iterations (def $NSmooth)"
  echo "   -noseghead : do not seg the head, just tess and smooth existing"
  echo "   -thresh1 threshold : default is $thresh1"
  echo "   -thresh2 threshold : default is $thresh2"
  echo "   -nhitsmin N        : default is $nhitsmin"
  echo "   -ndilate  N        : default is $ndilate"
  echo "   -nerode   N        : default is $nerode"
  echo "   -fillval value     : default is $fillval"
  echo "   -fhi fhi           : fhi for MRIchangeType(); default is to use default in mri_seghead"
  echo "   -no-rescale        : do not  rescale input when converting to uchar (-no-rescale)"
  echo "   -no-fill-holes-islands : do not fill holes and remove islands"
  echo "   -or-mask ormask.mgz : include all voxels in ormask in the head seg (-no-or-mask)"
  echo "   -tess/-mc : tessellation method using mri_tessellate or mri_mc (default is -mc)"
  echo "   -inflate : inflate and compute sulc"
  echo "   -curv : compute curv with smoothing"
  echo "   -srcvol volid      : default is T1"
  echo "   -headvol volid     : default is seghead"
  echo "   -headsurf surfid   : default is seghead"
  echo "   -smheadsurf surfid : default is smseghead"
  echo "   -hemi hemi         : default is lh"
  echo "   -sd subjectsdir    : default is SUBJECTS_DIR"
  echo "   -umask umask       : default is 2 (ie, group and individ writable)"
  echo "   -log logfile"
  echo ""
  echo "See also: mri_seghead, mri_tessellate, and mris_smooth"
  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 will segment and create a surface representation of the head that
can be displayed using tksurfer. The output is lh.seghead and will
be stored in the subject's surf directory. There will also be a new
volume called seghead in the subject's mri directory. The volume will
be segmented and filled with 255. If the final head surface does not
look good enough, this volume can be edited (eg, tkmedit subjid
seghead -aux T1), then mkheadsurf can be re-run with -noseghead. The
segmenation surface can be loaded by File->Load Main Surface:
lh.smseghead.

For notes on setting the head segmentation parameters thresh1,
thresh2, and nhitsmin, see mri_seghead --help.

tksurfer subject lh seghead