#! /bin/tcsh -f

#
# make_average_surface
#
# Creates average surfaces and curvatures from a set of subjects.
#
# --help option will show usage
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: greve $
#    $Date: 2009/09/21 17:14:26 $
#    $Revision: 1.51 $
#
# Copyright (C) 2002-2010,
# 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
#

# uncomment this to increase number of allowable open files to maximum:
#limit descriptors unlimited

set VERSION = '$Id: make_average_surface,v 1.51 2009/09/21 17:14:26 greve Exp $';
setenv FIX_VERTEX_AREA

set PrintHelp = 0;
set UseSymLink = 1;
set sdout = ();
set inputargs = ($argv);

if ( $?SUBJECTS) then
  set SUBJECTS=($SUBJECTS)
endif
if ( $?SUBJECTS_DIR) then
  set SUBJECTS_DIR=($SUBJECTS_DIR)
endif
set average_subject=average
set sphere_reg=sphere.reg
set transform_fname=""
set XFORM_FLAG=""
set ico_no=7
set DoXHemi = 0;
set topdir = ();
set inputargs = ($argv)

# these should match latest declarations in recon-all
set DESIKAN_GCS = curvature.buckner40.filled.desikan_killiany.2010-03-25.gcs
set DESTRIEUX_GCS = destrieux.simple.2009-07-29.gcs

# note: SUBJECTS can be set by --subjects
# SUBJECTS_DIR can be set by --sdir
# average_subject can be set by --out
# sphere_reg can be set by --surf-reg
# transform_fname can be set by --xform (default is no file specified)
# XFORM_FLAG appends -X to transform_fname if --xform is specified

if($#argv == 0) then
  # zero args is allowed only if SUBJECTS env var is declared
  if ( ! $?SUBJECTS) then
    goto usage_exit;
  endif
endif

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:


# This will put the data under $sdout (SUBJECTS_DIR by default)
set ddir = ${sdout}/${average_subject}

mkdir -p ${ddir} \
        ${ddir}/surf \
        ${ddir}/mri \
        ${ddir}/mri/transforms \
        ${ddir}/scripts \
        ${ddir}/label \
        ${ddir}/stats

set LF = ${ddir}/scripts/make_average_surface.log
if(-e $LF) mv $LF $LF.bak

echo "" | tee -a $LF
echo $0 | tee -a $LF
echo $inputargs | tee -a $LF
echo "" | tee -a $LF
echo $VERSION | tee -a $LF
date | tee -a $LF
pwd | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "output ddir is $sdout" | tee -a $LF
id | tee -a $LF

#
# Begin script guts...
#
echo ==================== | tee -a $LF
echo make_average_surface | tee -a $LF
echo ==================== | tee -a $LF
date | tee -a $LF
pwd | tee -a $LF
echo $0 | tee -a $LF
echo $inputargs | tee -a $LF
echo input subjects: ${SUBJECTS}  | tee -a $LF
echo output subject: ${average_subject}  | tee -a $LF

# Copy MNI 305 to output mri dir
set mni305 = $FREESURFER_HOME/average/mni305.cor.mgz
if(! -e $mni305) then
  echo "ERROR: cannot find $mni305"
  exit 1;
endif
cp $mni305 $sdout/${average_subject}/mri

# Set the MNI305 xfm fname to 'auto'
set mnivol = $sdout/${average_subject}/mri/mni305.cor.mgz
set cmd = (mri_add_xform_to_header -c auto $mnivol $mnivol)
echo $cmd | tee -a $LF
$cmd  | tee -a $LF
if($status) exit 1;

# Include cross-hemi registration
if($DoXHemi) then
  set tmp = ();
  foreach s ($SUBJECTS)
    set tmp = ($tmp $s $s/xhemi)
  end
  set SUBJECTS = ($tmp);
endif

set tmp_subjects = ( $SUBJECTS )

#------------------------------------------------------------
foreach hemi (lh rh)

  echo ""  | tee -a $LF
  echo ""  | tee -a $LF

  cd $sdout/${average_subject}

  # ------------------------------------------------------------------
  echo "#@# Making $hemi registration template ---------"  | tee -a $LF
  set atlas = $hemi.reg.template.tif
  set cmd = (mris_make_template -norot -annot aparc \
             $hemi $sphere_reg $tmp_subjects $atlas)
  pwd | tee -a $LF
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit(1);

  #--------------------------------------------------------------
  # Do we really want to do this?
  echo "#@# Making $hemi annotation ---------------------"  | tee -a $LF
  date  | tee -a $LF
  set GCS     = $DESIKAN_GCS
  set CPAtlas = $hemi.$GCS
  set annot = label/$hemi.aparc.annot
  set cmd = (mris_ca_label -ml-annot $CPAtlas $ico_no $annot)
  echo $cmd  | tee -a $LF
  pwd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) then
    echo "ERROR: mris_ca_label failed"  | tee -a $LF
    exit 1;
  endif

  set GCS     = $DESTRIEUX_GCS
  set CPAtlas = $hemi.$GCS
  set annot   = label/$hemi.aparc.a2009s.annot
  set cmd = (mris_ca_label -ml-annot $CPAtlas $ico_no $annot)
  echo $cmd  | tee -a $LF
  pwd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) then
    echo "ERROR: mris_ca_label failed"  | tee -a $LF
    exit 1;
  endif

  #--------------------------------------------------------------
  foreach s (orig white pial inflated)
    echo "#@# Making average $hemi.$s surface ---------------------"  | tee -a $LF
    date  | tee -a $LF
    set cmd = (mris_make_average_surface $XFORM_FLAG \
                -i $ico_no -o $s -sdir-out $sdout\
                ${hemi} \
                ${s}_avg ${sphere_reg} \
                ${average_subject} \
                $SUBJECTS )
    pwd  | tee -a $LF
    echo $cmd  | tee -a $LF
    $cmd |& tee -a $LF
    if($status) then
      echo "mris_average_surface failed"  | tee -a $LF
      exit 1;
    endif
    echo ""  | tee -a $LF
    echo ""  | tee -a $LF
  end
  # mris_make_average_surface creates its own sphere.reg, so copy over it with ico
  mv surf/${hemi}.${sphere_reg} surf/${hemi}.${sphere_reg}.avg
  if($ico_no == 7) then
    # for historical reasons, use ?h.sphere.reg for ico7
    cp ${FREESURFER_HOME}/average/surf/${hemi}.sphere.reg surf/${hemi}.${sphere_reg}
    cp ${FREESURFER_HOME}/average/surf/${hemi}.sphere.reg surf/${hemi}.sphere.reg
  else
    cp ${FREESURFER_HOME}/average/surf/${hemi}.sphere.ico$ico_no.reg surf/${hemi}.${sphere_reg}
    cp ${FREESURFER_HOME}/average/surf/${hemi}.sphere.ico$ico_no.reg surf/${hemi}.sphere.reg
  endif
  echo "---------------------------------------------------"  | tee -a $LF
  echo "" | tee -a $LF
  echo "" | tee -a $LF
  echo "" | tee -a $LF

  #--------------------------------------------------------------
  foreach c (sulc curv thickness)
    echo "#@# Making average $hemi.$c ----------------------" | tee -a $LF
    date
    set cmd = (mris_average_curvature \
                -o ${average_subject} \
                $c \
                ${hemi} \
                ${sphere_reg} \
                $SUBJECTS \
                $sdout/${average_subject}/surf/${hemi}.avg_${c})
    pwd | tee -a $LF
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) then
      echo "mris_average_curvature failed" | tee -a $LF
      exit 1;
    endif
    echo "" | tee -a $LF
    echo "" | tee -a $LF

    if(0 && $UseSymLink) then
      pushd $sdout/$average_subject/surf/
      ln -s -f $hemi.avg_$c $hemi.$c
      popd
    else
      cp $sdout/$average_subject/surf/$hemi.avg_$c \
         $sdout/$average_subject/surf/$hemi.$c
    endif

  end
  echo "---------------------------------------------------" | tee -a $LF
  echo "" | tee -a $LF
  echo "" | tee -a $LF
  echo ""

  #--------------------------------------------------------------
  echo "#@# Smoothing, inflating, etc $hemi ----------------" | tee -a $LF
  cd $sdout/${average_subject}/surf
  pwd     | tee -a $LF

  # This writes out ?h.area and ?h.curv. Want area, not curv.
  set cmd = (mris_smooth -n 5 ./${hemi}.white_avg ./${hemi}.smoothwm)
  echo $cmd |& tee -a $LF
  $cmd | tee -a $LF
  echo "" | tee -a $LF
  rm $hemi.curv
  cp $hemi.avg_curv $hemi.curv

  set niters = "";
  if($ico_no <= 5) set niters = "-n 2";
  set cmd = (mris_inflate -no-save-sulc $niters ./${hemi}.smoothwm ./${hemi}.inflated)
  echo $cmd |& tee -a $LF
  $cmd | tee -a $LF
  echo "" | tee -a $LF

  set cmd = (mris_curvature -a 10 -w ./${hemi}.white_avg)
  echo $cmd |& tee -a $LF
  $cmd | tee -a $LF
  echo "" | tee -a $LF

  if($UseSymLink) then
    ln -s -f ${hemi}.${sphere_reg} ${hemi}.sphere
    ln -s -f ${hemi}.orig_avg  ${hemi}.orig
    ln -s -f ${hemi}.white_avg ${hemi}.white
    ln -s -f ${hemi}.pial_avg  ${hemi}.pial
  else
    cp -f ${hemi}.${sphere_reg} ${hemi}.sphere
    cp -f ${hemi}.orig_avg  ${hemi}.orig
    cp -f ${hemi}.white_avg ${hemi}.white
    cp -f ${hemi}.pial_avg  ${hemi}.pial
  endif

  # -----------------------------------------------------
  # Create cortex label from Christophe's medial wall
  set annot = $ddir/label/$hemi.aparc.a2009s.annot
  if(! -e $annot) then
    echo "ERROR: cannto find $annot" |& tee -a $LF
    exit 1;
  endif
  # Create tmpdir
  set tmpdir = $ddir/make_cortex_label.tmp.$$
  rm -rf $tmpdir
  mkdir -p $tmpdir
  # Unload all of Christophe's labels into tmpdir
  set cmd = (mri_annotation2label --subject $average_subject \
    --hemi $hemi --outdir $tmpdir --annotation aparc.a2009s)
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  mv $tmpdir/$hemi.Medial_wall.label $ddir/label; # Move medial wall
  # Create the cortex label with all the remaining labels
  set ctxlabel = $ddir/label/$hemi.cortex.label
  set cmd = (mri_mergelabels -o $ctxlabel -d $tmpdir)
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  rm -rf $tmpdir; # Cleanup
  # -------- done creating cortex label -----------------

end

#--------------------------------------------------------------
# Now create the average area surface
set subjlistfile = /tmp/subjlist.$$
rm -f $subjlistfile
echo $SUBJECTS > $subjlistfile
foreach hemi (rh lh)
  foreach surf (orig white pial)
    set avgarea = $hemi.$surf.avg.area.mgh
    set cmd = (mris_preproc --out $avgarea  --f $subjlistfile)
    set cmd = ($cmd --target $average_subject --hemi $hemi)
    set cmd = ($cmd --area $surf --mean --surfreg $sphere_reg)
    pwd | tee -a $LF
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit (1)
  end
end

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

exit 0

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

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

  set flag = $argv[1];
  if (! $getting_subjects) then
    shift;
  endif

  switch($flag)

    case "--help":
      set PrintHelp = 1;
      goto usage_exit;
      exit 0;

    case "--version":
      echo $VERSION
      exit 0;

    case "--subjects":
      if ( $#argv == 0) goto arg1moreerr;
      set SUBJECTS = $argv[1]; shift;
      # loop on getting variable number of subject names
      set getting_subjects = 1; # see 'default:' case
      breaksw

    case "--fsgd":
      if ( $#argv == 0) goto arg1err;
      set fsgdf = $argv[1]; shift;
      if(! -e $fsgdf) then
        echo "ERROR: cannot find $fsgdf";
        exit 1;
      endif
      set sl = `cat $fsgdf | awk '{if($1 == "Input") print $2}'`;
      set SUBJECTS = ($sl);
      breaksw

    case "--out":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if ( $#argv == 0) goto arg1err;
      set average_subject = $argv[1]; shift;
      breaksw

    case "--surfreg":
    case "--surf-reg":
    case "--surf_reg":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if ( $#argv == 0) goto arg1err;
      set sphere_reg = $argv[1]; shift;
      breaksw

    case "--topdir":
    case "--sd-out":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if( $#argv < 1) goto arg1err;
      set sdout = $argv[1]; shift;
      breaksw

    case "--sd":
    case "--sdir":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if ( $#argv == 0) goto arg1err;
      set SUBJECTS_DIR = $argv[1]; shift;
      setenv SUBJECTS_DIR $SUBJECTS_DIR
      breaksw

    case "--xform":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if ( $#argv == 0) goto arg1err;
      set transform_fname = $argv[1]; shift;
      set XFORM_FLAG = "-X $transform_fname";
      breaksw

    case "--ico":
    case "--i":
      if ( $getting_subjects ) then
        # got em all, from --subjects variable arg loop
        set getting_subjects = 0;
        shift;
      endif
      if ( $#argv == 0) goto arg1err;
      set ico_no = $argv[1]; shift;
      breaksw

    case "--symlink":
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      set UseSymLink = 1;
      breaksw

    case "--no-symlink":
    case "--no-link":
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      set UseSymLink = 0;
      breaksw

    case "--nocleanup":
      # This has no effect, but needed for compat with make_average_volume
      set cleanup = 0;
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      breaksw

    case "--debug":
    case "--echo":
      set echo = 1;
      set verbose = 1;
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      breaksw

    case "--xhemi":
      set DoXHemi = 1;
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      breaksw

    case "--no-xhemi":
      set DoXHemi = 0;
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      breaksw

    # These are flags passed to make_average_subject, but dont apply here
    case "--no-surf":
    case "--no-vol":
    case "--force":
    case "--link":
    case "--keep-all-orig":
      if ( $getting_subjects ) then
        set getting_subjects = 0;
        # got em all, from --subjects variable arg loop
      endif
      breaksw

    default:
      if ( $getting_subjects ) then
        # loop on getting variable number of subjects,
        # until a new flag is found, or no more args
        set SUBJECTS = ( $SUBJECTS $argv[1] ); shift;
        set getting_subjects = 1;
      else
        echo ERROR: Flag $flag unrecognized.
        echo $cmdline
        exit 1
      endif
      breaksw
  endsw

end

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

############--------------##################
check_params:
  if (! $?SUBJECTS) then
    echo "ERROR: no subjects declared!"
    echo "Either declare subjects in SUBJECTS variable,"
    echo "or declare using --subjects argument."
    exit 1
  endif
  if (! $?SUBJECTS_DIR) then
    echo "ERROR: SUBJECTS_DIR is not declared!"
    echo "Either set the SUBJECTS_DIR environment variable,"
    echo "or declare using --sdir argument, the root directory"
    echo "for subject data files."
    exit 1
  endif
  if(! -e $SUBJECTS_DIR ) then
    echo "ERROR: SUBJECTS_DIR $SUBJECTS_DIR does not exist."
    exit 1;
  endif
  if(! $?FREESURFER_HOME ) then
    echo "ERROR: environment variable FREESURFER_HOME not set."
    exit 1;
  endif
  if(! -e $FREESURFER_HOME ) then
    echo "ERROR: FREESURFER_HOME $FREESURFER_HOME does not exist."
    exit 1;
  endif
  if($#sdout == 0) set sdout = $SUBJECTS_DIR;
goto check_params_return;
############--------------##################

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

############--------------##################
arg1moreerr:
  echo "ERROR: flag $flag requires one or more arguments"
  exit 1
############--------------##################

############--------------##################
usage_exit:
  echo ""
  echo "USAGE: make_average_surface"
  echo ""
  echo "Required Arguments"
  echo "   --subjects <subj1> <subj2> ... <subjN>"
  echo "             : or declare subjects in SUBJECTS env var"
  echo "   --fsgd fsgdfile : get subject list from fsgd"
  echo ""
  echo "Optional Arguments"
  echo "   --out <average subject name>    : default name is 'average'"
  echo "   --sdir <SUBJECTS_DIR to use instead of the one in the env>"
  echo "   --sd      : same as --sdir"
  echo "   --sd-out topdir : put data here"
  echo "   --xform <transform_fname>       : filename of transform file"
  echo "   --ico <number>  : specify icosahedron number, default is 7 (ic7.tri)"
  echo "   --surf-reg <surface name>   : alternative registration surface"
  echo "                                 default: sphere.reg"
  echo "   --force   : overwrite existing average subject data"
  echo ""
  echo "   --help       : short descriptive help"
  echo "   --no-symlink : do not use symbolic links (just copy files)"
  echo "   --no-lik     : same as --no-symlink"
  echo "   --version : script version info"
  echo "   --echo    : enable command echo, for debug"
  echo "   --debug   : same as --echo"
  echo ""
  echo "See also: recon-all, make_final_surfaces, morph_subject"
  echo ""

  if(! $PrintHelp) exit 1;

  echo Version: $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

Creates average surfaces and curvatures from a set of subjects.

Calls mris_average_curvature, mris_make_average_surface, mris_smooth,
mris_inflate, and mris_curvature.

EXAMPLE

make_average_surface --out avgsubject --subjects subj1 subj2 subj3 subj4

will create $SUBJECTS_DIR/avgsubject with average surfaces for orig,
white, pial, inflated for each hemi.

SEE ALSO

make_average_subject, make_average_volume

GETTING HELP

Run recon-all --help for extensive text on the reconstruction process.
Or, send email to freesurfer@nmr.mgh.harvard.edu.
Also see http://surfer.nmr.mgh.harvard.edu.

