#!/bin/tcsh -f

@global_parse `basename $0` "$*" ; if ($status) exit 0

# make surface spec files from the surface files
#
# usage @SUMA_Make_Spec_FS [options]
#
# options:  -sid   subject_id   : specify subject ID
#       -fspath   FreeSurfer_path : specify path to FreeSurfer files
#       -neuro         : use neurological orientation
#           -debug      level      : display extra output
#

#----------------------------------------------------------------------
goto L_INIT_VARS
     L_INIT_VARS_DONE:

goto L_CHECK_USAGE
     L_CHECK_USAGE_DONE:

goto L_PARSE_COMMAND
     L_PARSE_COMMAND_DONE:

goto L_VERIFY_PROGRAMS
     L_VERIFY_PROGRAMS_DONE:

goto L_SET_SURF_DIRS
     L_SET_SURF_DIRS_DONE:

goto L_CHECK_FOR_OVERWRITE
     L_CHECK_FOR_OVERWRITE_DONE:

goto L_LOOK_FOR_SURF
     L_LOOK_FOR_SURF_DONE:

goto L_CREATE_BRICK
     L_CREATE_BRICK_DONE:

goto L_CONVERT_SURFACES
     L_CONVERT_SURFACES_DONE:

goto L_CREATE_SPEC
     L_CREATE_SPEC_DONE:

goto L_ICO
     L_ICO_DONE:

goto L_TEST_SURF_VOL
     L_TEST_SURF_VOL_DONE:

goto L_GOOD_END   # finished, woohooo!


####################################################################
# variable initialization

L_INIT_VARS:

    set prog_name = $0:t
    set endstr = "$prog_name ... finished"
    set debug     = 0
    set ldlist = (141 60)
    #set fsdsets = (thickness.gii.dset)
    set fsdsets = (thickness curv sulc) 
    set inflates = ()
    
    set surf_attribs = (smoothwm      \
         pial         \
         inflated      \
         occip.patch.3d      \
         occip.patch.flat   \
         occip.flat.patch.3d \
         fusiform.patch.flat  \
         full.patch.3d      \
         full.patch.flat   \
         full.flat.patch.3d \
         full.flat   \
         flat.patch      \
         sphere         \
         white          \
         sphere.reg     \
         rh.sphere.reg  \
         lh.sphere.reg  \
         pial-outer-smoothed     \
         )

goto L_INIT_VARS_DONE


####################################################################
# check usage, and possibly print help

L_CHECK_USAGE:

    if ( $#argv == 0 ) then
   echo "usage: $prog_name [options] -sid SUBJECT_ID"
   echo "usage: $prog_name -help"
   set endstr = ""
   goto L_GOOD_END
    endif

goto L_CHECK_USAGE_DONE


####################################################################
# parse the command line

L_PARSE_COMMAND:

    # init command line arg values

    set fs_dir      = "."
    set subj_id      = ""
    set neuro_ori   = 0
    set sfieldname  = 'FreeSurferSurface'
    set sex = 'asc'
    set stp = 'FreeSurfer'
    set ldu = ()
    set ldpref = ()
    set use_mgz = 0
    set set_space = ''
    set args     = $#argv
    set count    = 1

    while ( $count <= $args )
   switch ( "$argv[$count]" )

       # ----------------------------------------------------------
       # usage: -help
       case "-h":
       case "-help":

      goto L_HELP_END      # and don't ya' come back, neither

       breaksw

       # ----------------------------------------------------------
       # usage: -sid SUBJECT_ID
       case "-sid":

      if ( $count >= $args ) then
          set endstr = "arg usage: -sid SUBJECT_ID"
          goto L_BAD_END
      endif

      @ count ++
      set subj_id  = $argv[$count]

       breaksw

       # ----------------------------------------------------------
       # usage: -fspath FREESURFER_PATH
       case "-fspath":

      if ( $count >= $args ) then
          set endstr = "arg usage: -fspath FREESURFER_PATH"
          goto L_BAD_END
      endif

      @ count ++
      set fs_dir  = $argv[$count]

      if ( ! -d $fs_dir ) then
          set endstr = "failure: directory not found - '$fs_dir'"
          goto L_BAD_END
      endif

       breaksw

       # ----------------------------------------------------------
       # usage : -neuro
       case "-neuro":

      set neuro_ori = 1
       breaksw

       # ----------------------------------------------------------
       # usage : -nocor
       case "-nocor":
      set endstr = "Option -nocor obsolete. See -help for important details."
      goto L_BAD_END
      
      #set neuro_ori = -1
      #set sfieldname = 'SurfaceName'
      #set sex = 'gii'
      #set stp = 'GIFTI'
      
       breaksw

       # ----------------------------------------------------------
       # usage : -GNIFTI
       case "-GNIFTI":
       case "-NIFTI":
       case "-GIFTI":
       case "-IFTI":

      set neuro_ori = -1
      set sfieldname = 'SurfaceName'
      set sex = 'gii'
      set stp = 'GIFTI'
      
       breaksw

       # ----------------------------------------------------------
       # usage : -set_space SPACE
       case "-set_space":

      if ( $count >= $args ) then
          set endstr = "arg usage: -set_space SPACE"
          goto L_BAD_END
      endif

      @ count ++
      set set_space  = $argv[$count]

       breaksw

       # ----------------------------------------------------------
       # usage : -neuro
       case "-use_mgz":

      set use_mgz = 1

       breaksw

       # ----------------------------------------------------------
       # usage : -debug DEBUG_LEVEL
       case "-debug":

      if ( $count >= $args ) then
          set endstr = "arg usage: -debug DEBUG_LEVEL"
          goto L_BAD_END
      endif

      @ count ++
      set debug = $argv[$count]

      if ( "$debug" > 2 ) set debug = 2
      if ( "$debug" < 0 ) set debug = 0

       breaksw

       # ----------------------------------------------------------
       # usage : -inflate VAL
       case "-inflate":

      if ( $count >= $args ) then
          set endstr = "arg usage: -inflate INF"
          goto L_BAD_END
      endif

      @ count ++
       set ff = 0
       foreach inf ($inflates)
         if ($inf == "$argv[$count]") set ff = 1
       end
       if ($ff == 0) set inflates = ($inflates $argv[$count])

       breaksw


       # ----------------------------------------------------------
       # usage : -ld VAL
       case "-ld":

      if ( $count >= $args ) then
          set endstr = "arg usage: -ld LD"
          goto L_BAD_END
      endif

      @ count ++
      set ldu = ($ldu $argv[$count])

       breaksw

       # ----------------------------------------------------------
       # usage : -ldpref LDpref
       case "-ldpref":

      if ( $count >= $args ) then
          set endstr = "arg usage: -ldpref LDpref"
          goto L_BAD_END
      endif

      @ count ++
      set ldpref = ($ldpref $argv[$count])

       breaksw

       # ----------------------------------------------------------
       # usage : -no_ld
       case "-no_ld":

      set ldlist = ()

       breaksw

      
       # ----------------------------------------------------------
       # bad argument
       default:

      set endstr = "invalid option: '$argv[$count]'"
      apsearch -popt `basename $0` -word $argv[$count]
      goto L_BAD_END

       breaksw
   endsw

   @ count ++
    end

    if ( $subj_id == "" ) then
   set endstr = "missing required option: -sid"
   goto L_BAD_END
    endif

    if ($#inflates == 0 && $neuro_ori == -1) set inflates = ($inflates 200) 

    if ( $debug ) echo "-- usage okay"
    if ( $debug > 1 ) set echo
    if ( $#ldu ) set ldlist = ($ldu)
    if ($#ldpref == 0 && $#ldlist != 0) then
      set ldpref = ()
      foreach lldd ($ldlist)
         set ldpref = ($ldpref std.$lldd.)
      end
    endif
    if ($#ldpref != 0 && $#ldpref != $#ldlist) then
      set endstr = "Not enough -ldpref options for $#ldlist values"  
      goto L_BAD_END
    endif

    set spec_files   = ({$subj_id}_lh.spec {$subj_id}_rh.spec )
    set afni_prefix  = ${subj_id}_SurfVol
      if ( $neuro_ori != -1 ) then
    set afni_dataset = $afni_prefix+orig
      else
    set afni_dataset = $afni_prefix 
      endif
    # set and go to the base directory
    cd $fs_dir
    set start_dir = $cwd

    if ( $debug ) echo "-- using start_dir '$start_dir'"

goto L_PARSE_COMMAND_DONE


####################################################################
# make sure programs exist

L_VERIFY_PROGRAMS:

    set failed_pgms = 0
    foreach test_command ( afni to3d suma mris_convert )

   (which $test_command) >& /dev/null

   if ( $status ) then
       echo "program not found in path: $test_command"
       @ failed_pgms ++
   endif
    end

    if ( $failed_pgms ) then
   set endstr = "$failed_pgms program(s) not found"
   goto L_BAD_END
    endif

    if ( $debug ) echo "-- all programs found"

goto L_VERIFY_PROGRAMS_DONE


####################################################################
# 1. set surf_dir and orig_dir - check cwd and parent
# 2. create SUMA directory at the same level as surf_dir and store the results there

L_SET_SURF_DIRS:

    # find surf directory

    if ( -d surf ) then
   set surf_dir = ./surf
   set label_dir = ./label
    else if ( $cwd:t == surf ) then
   set surf_dir = .
   set label_dir = ../label
    else if ( -d ../surf ) then
   set surf_dir = ../surf
   set label_dir = ../label
    else
   # this is a general failure case, even if we find one
   set surf_dirs = ( `find . -maxdepth 4 -type d -name surf -print` )
   
   if ( $#surf_dirs == 0 ) then
       echo "failure: cannot find directory 'surf' under '$fs_dir'"
       echo "(subject to a maximum search depth of 4 subdirectories)"

       set endstr = ""
   else if ( $#surf_dirs == 1 ) then
       echo "surf directory found at '$surf_dirs[1]'"
       set endstr = "consider running program from '$surf_dirs[1]:h'"
   else
       echo "multiple surf directories found:"
       set count = 1
       while ( $count <= $#surf_dirs )
      echo "     $surf_dirs[$count]"
      @ count ++
       end

       set endstr = ( "consider running program from one of the" \
            "surf directories" )
   endif

   goto L_BAD_END
    endif

    # verify surf dir permissions
    if ( ! -w $surf_dir ) then
   set endstr = "failure: no write permissions for directory '$surf_dir'"
   goto L_BAD_END
    endif

    if ( $debug ) echo "-- using surf directory '$surf_dir'..."

    # now check for orig dir

    set orig_dir = ""
    foreach test_dir (  $surf_dir/orig      \
         $surf_dir/mri/orig   \
         orig         \
         mri/orig      \
         ../orig         \
         ../mri/orig )
   if ( -d $test_dir ) then
       set orig_dir = $test_dir
       break
   endif
    end

    set orig_mgz = ""
    foreach test_mgz (  $surf_dir/orig.mgz      \
         $surf_dir/mri/orig.mgz   \
         orig.mgz         \
         mri/orig.mgz      \
         ../orig.mgz         \
         ../mri/orig.mgz )
   if ( -f $test_mgz ) then
       set orig_mgz = $test_mgz
       break
   endif
    end
    
   set other_mgz = ()
   set other_candidates = ( T1 aparc+aseg aparc.a2005s+aseg aparc.a2009s+aseg aseg.auto aseg brain.finalsurfs brain brainmask.auto brainmask norm nu nu_noneck lh.ribbon rh.ribbon wm.asegedit wm wm.seg )
   foreach candidate ($other_candidates)
      foreach dcand ($surf_dir/ $surf_dir/mri/ ./ mri/ ../ ../mri/ )
         if ( -f $dcand$candidate.mgz ) then 
            set other_mgz = ($other_mgz $dcand$candidate.mgz)
            break
         endif
      end
   end
   
    if ( $orig_dir == "" && $orig_mgz == "") then
   set endstr = "failure: cannot find directory 'orig' or file 'orig.mgz' under $fs_dir"
   goto L_BAD_END
    endif

    if ( $orig_dir == "" ) then
   set endstr = "failure: script expects an 'orig' directory even if it is empty."
   goto L_BAD_END 
    endif
    
    # verify orig dir permissions
    if ( ! -w $orig_dir ) then
   set endstr = "failure: no write permissions for directory '$orig_dir'"
   goto L_BAD_END
    endif

    if ( $debug ) echo "-- using orig directory '$orig_dir'..."

    # decide whether we need to make surf/SUMA directory

    set suma_dir = ./SUMA

    if ( -d $suma_dir ) then
   if ( $debug ) echo "-- $suma_dir already exists, continuing..."
    else
   echo "++ creating directory '$suma_dir' for results..."
   mkdir $suma_dir
   if ( $status ) then
       set endstr = "failure: cannot create directory '$suma_dir'"
       goto L_BAD_END
   endif
    endif

goto L_SET_SURF_DIRS_DONE


####################################################################
# verify non-existence of spec files and AFNI files

L_CHECK_FOR_OVERWRITE:

    set test_failures = 0

    foreach test_file (   $suma_dir/$spec_files[1]   \
         $suma_dir/$spec_files[2]   \
         $suma_dir/$afni_dataset.HEAD   \
         $suma_dir/$afni_dataset.BRIK   \
         $orig_dir/$afni_dataset.HEAD   \
         $orig_dir/$afni_dataset.BRIK   \
         $orig_dir/$afni_prefix'+orig.HEAD' \
         $orig_dir/$afni_prefix'+orig.BRIK' \
         $suma_dir/$afni_prefix'+orig.HEAD' \
         $suma_dir/$afni_prefix'+orig.BRIK' \
         $suma_dir/$afni_prefix'.nii' \
         $suma_dir/$afni_prefix'.nii.gz')
   if ( -f $test_file ) then
       if ( $test_failures == 0 ) then
      echo "failure: will not overwrite files: "
      set test_failures = 1
       endif

       echo "        '$test_file'"
   endif
    end

    if ( $test_failures ) then
        set endstr = "please remove these files if you want to rerun the script."
   goto L_BAD_END
    endif

goto L_CHECK_FOR_OVERWRITE_DONE

####################################################################
# find surface files

L_LOOK_FOR_SURF:

    cd $surf_dir

    set list_lh = "X"   # init to something useless - allows nice empty check
    set list_rh = "X"

    foreach attrib ( $surf_attribs )
   if ( -f lh.$attrib ) set list_lh = ( $list_lh lh.$attrib )
   if ( -f rh.$attrib ) set list_rh = ( $list_rh rh.$attrib )
    end

    set list_lh = ( $list_lh[2-] )   # now remove the leading "X"
    set list_rh = ( $list_rh[2-] )

    if ( $#list_lh == 0 && $#list_rh == 0 ) then
   set endstr = "found no LH or RH surface files under '$surf_dir'"
   goto L_BAD_END
    endif

    if ( $#list_lh > 0 ) then
   echo "-- found $#list_lh LH surfaces"
   if ( $debug ) echo "    --" $list_lh
    endif

    if ( $#list_rh > 0 ) then
   echo "-- found $#list_rh RH surfaces"
   if ( $debug ) echo "    --" $list_rh
    endif

    cd $start_dir

goto L_LOOK_FOR_SURF_DONE


####################################################################
# convert to ascii via mris_convert (found earlier)

L_CONVERT_SURFACES:

    cd $surf_dir
    foreach surf ( $list_lh $list_rh )
   # run mris_convert and verify
   if ( $surf =~ *patch* || $surf =~ *.flat) then
       echo  "-- running 'mris_convert -p $surf $surf.asc'"
       mris_convert -p $surf $surf.asc
   else
       echo "-- running 'mris_convert $surf $surf.asc'"
       mris_convert $surf $surf.asc >& ___out
       tail ___out ; \rm -f ___out
   endif

   if ( $neuro_ori == -1 ) then
      if ( ! -f xmatras.1D ) then
         set c_ras = (`mri_info --cras --ras_good $start_dir/$orig_mgz`)
         set minus = `ccalc "$c_ras[3]*-1"` 
         echo 1 0 0 $c_ras[1] 0 1 0 $c_ras[2] 0 0 1 $minus > xmatras.1D
      endif
      #If using new conversion, change surfaces to gifti with proper
      #coordinate shift.
      if ( $surf =~ *.white* ||     \
           $surf =~ *.smoothwm* ||  \
           $surf =~ *.pial*) then
         ConvertSurface -sv $start_dir/$suma_dir/$afni_dataset.nii -i $surf.asc \
                        -ixmat_1D xmatras.1D -overwrite -o $surf.gii
      else
         ConvertSurface -sv $start_dir/$suma_dir/$afni_dataset.nii -i $surf.asc \
                        -ixmat_1D xmatras.1D -overwrite -o $surf.gii
      endif
   endif

   if ( ! -f $surf.${sex} ) then
       echo "failure: could not create surface $surf.${sex}"

       if ( ! -w . ) then
      set endstr = "--> no write permissions in this directory"
       else
      set endstr = "--> is your FreeSurfer license installed?"
       endif

       goto L_BAD_END
   endif

   if ( $debug ) then
       echo "++ $surf.${sex} created"
       if ( -f $start_dir/$suma_dir/$surf.${sex} ) then
      echo "++ overwriting '$suma_dir/$surf.${sex}'"
       endif
   endif
   
   #tmp fix for rh.lh.sphere.reg.gii side 
   if ( "$surf.${sex}" == rh.lh.sphere.reg.gii) then
      echo "++ Fixing side for rh.lh.sphere.reg.gii"
      sed 's:CortexLeft:CortexRight:' rh.lh.sphere.reg.gii > __tmpfix
      if ($status) then
         echo "** Failed to fix side for rh.lh.sphere.reg.gii"
      else 
         mv __tmpfix rh.lh.sphere.reg.gii
      endif 
   endif
   
   \mv -f $surf.${sex} $start_dir/$suma_dir   # how barbaric!

   if ( $status ) then
       set endstr = ( "failure: cannot write"   \
            "'$suma_dir/$surf.${sex}'" )
       goto L_BAD_END
   endif
    end

    cd $start_dir

   #and do the labels
   if ( -d $label_dir) then
      cd $label_dir
      foreach hand ( lh rh )
        if ( -f ${hand}.aparc.a2005s.annot ) then
          FSread_annot   -input ${hand}.aparc.a2005s.annot \
            -roi_1D  $start_dir/$suma_dir/${hand}.aparc.a2005s.annot.1D.roi \
            -dset    $start_dir/$suma_dir/${hand}.aparc.a2005s.annot.niml.dset \
            -cmap_1D $start_dir/$suma_dir/${hand}.aparc.a2005s.annot.1D.cmap 
        endif
     
        if ( -f ${hand}.aparc.a2009s.annot ) then
          FSread_annot   -input ${hand}.aparc.a2009s.annot \
            -roi_1D  $start_dir/$suma_dir/${hand}.aparc.a2009s.annot.1D.roi \
            -dset    $start_dir/$suma_dir/${hand}.aparc.a2009s.annot.niml.dset \
            -cmap_1D $start_dir/$suma_dir/${hand}.aparc.a2009s.annot.1D.cmap 
        endif

        foreach pp (prob predict)
           if ( -f ${hand}.v1.${pp}.label ) then
              sed '1,2s/^/#/' ${hand}.v1.${pp}.label > ${hand}.v1.${pp}.label.1D
              ConvertDset    -input ${hand}.v1.${pp}.label.1D'[$]' \
                             -node_index_1D ${hand}.v1.${pp}.label.1D'[0]' \
                             -o_niml \
                             -prefix $start_dir/$suma_dir/${hand}.v1.${pp} 
              \rm -f ${hand}.v1.${pp}.label.1D
              cd $start_dir/$suma_dir/
              3drefit -sublabel 0 v1_${pp} ${hand}.v1.${pp}.niml.dset
              cd -
           endif
        end
     
      end
      cd $start_dir
   endif
   
   #and the thickness data
   cd $surf_dir
   foreach hand ( lh rh )
      foreach tp ($fsdsets)
      if ( -f ${hand}.${tp} ) then
         mris_convert -c ${hand}.${tp} ${hand}.white ${hand}.${tp}.gii
         mv    ${hand}.${tp}.gii   \
               $start_dir/$suma_dir/${hand}.${tp}.gii.dset
      endif
      end
   end
   cd $start_dir
   
   #and make semi inflated surfaces
   cd $suma_dir
   if ($#inflates > 0) then
      foreach inf ($inflates)
         set surf_attribs = ($surf_attribs inf_${inf})
      end
      foreach hand ( lh rh )
         foreach inf ($inflates) 
            if ( -f ${hand}.white.asc ) then
               set iis = ${hand}.white.asc
            else if ( -f ${hand}.white.gii ) then
               set iis = ${hand}.white.gii
            else
               set endstr = ( "failure: Could not find white for inflation" )
               goto L_BAD_END
            endif
            SurfSmooth  -i ${iis}    -met NN_geom \
                        -Niter $inf  -o_gii   -surf_out ${hand}.inf_${inf} \
                        -match_vol 0.01
            if ( $hand == lh ) then
               set list_lh = ($list_lh ${hand}.inf_${inf})
            else
               set list_rh = ($list_rh ${hand}.inf_${inf})
            endif
         end
      end
   endif
   cd $start_dir
   
goto L_CONVERT_SURFACES_DONE


####################################################################
# actually create the spec file

L_CREATE_SPEC:

    cd $suma_dir

    foreach hand ( lh rh )
    if ( $hand == lh ) then
       set list_cur = ( $list_lh )      # get a current list copy
   else
       set list_cur = ( $list_rh )
   endif

   set spec_file = {$subj_id}_$hand.spec

   if ( $debug ) echo "++ creating spec file '$spec_file'..."

   (echo "# delimits comments" > $spec_file) >& /dev/null

   if ( $status ) then
       set endstr = ( "failure: no permissions to create spec file" \
          "'$suma_dir/$spec_file'" )
       goto L_BAD_END
   endif

   # note user, date, machine, pwd, command line
   echo ""                     >> $spec_file
   echo "# Creation information:"            >> $spec_file
   echo "#     user    : $user"            >> $spec_file
   echo "#     date    : `date`"            >> $spec_file
   echo "#     machine : `uname -n`"         >> $spec_file
   echo "#     pwd     : $cwd"            >> $spec_file
   echo "#     command : $prog_name $argv"         >> $spec_file
   echo ""                     >> $spec_file
   
   # define the group
   echo "# define the group"            >> $spec_file
   echo "        Group = $subj_id"            >> $spec_file
   echo ""                     >> $spec_file

   # define the states
   echo "# define various States"            >> $spec_file
   foreach attrib ( $surf_attribs )
       echo "        StateDef = $attrib"         >> $spec_file
   end
   echo ""                     >> $spec_file

   foreach surf ( $list_cur )
       set s_head = `echo $surf | cut -d. -f1`   # up to first '.'
       set s_state = `echo $surf | cut -d. -f2-`   # after first '.'

      set label_dset = ""      
       # check for SAME mapping ref
       if ( "$surf.${sex}" == "$s_head.smoothwm.${sex}" ) then
      set map_ref = SAME
         if ( -f ${hand}.aparc.a2005s.annot.niml.dset ) then
      set label_dset = ${hand}.aparc.a2005s.annot.niml.dset
         endif
         if ( -f ${hand}.aparc.a2009s.annot.niml.dset ) then
      set label_dset = ${hand}.aparc.a2009s.annot.niml.dset
         endif
       else
      set map_ref = $s_head.smoothwm.${sex}
       endif
       
       # check for SAME mapping ref case for pial-outer
       if ( "$surf.${sex}" == "$s_head.pial-outer-smoothed.${sex}" ) then
      set map_ref = SAME
       else
      #abide by previous setting
       endif
       
       # EmbedDimension is 2 for .flat surfaces, else 3
       if ( $surf =~ *.flat* ) then
      set embed_ref = 2
       else
      set embed_ref = 3
       endif

       # Anatomical flag
       if ( $surf =~ *.white* ||\
            $surf =~ *.smoothwm* ||\
            $surf =~ *.pial*) then
      set anatomical = Y
       else
      set anatomical = N
       endif
       
       echo "NewSurface"               >> $spec_file
       echo "        SurfaceFormat = ASCII"      >> $spec_file
       echo "        SurfaceType = ${stp}"      >> $spec_file
       echo "        ${sfieldname} = $surf.${sex}"   >> $spec_file
       echo "        LocalDomainParent = $map_ref"      >> $spec_file
       echo "        SurfaceState = $s_state"      >> $spec_file
       echo "        EmbedDimension = $embed_ref"      >> $spec_file
       echo "        Anatomical = $anatomical" >> $spec_file
         if ("$label_dset" != "") \
       echo "        LabelDset = $label_dset" >> $spec_file
       echo ""                  >> $spec_file
   end

   echo "++ created spec file'$suma_dir/$spec_file'"
    end    # foreach hand

   # Create both.spec file
   set spec_file = {$subj_id}_both.spec

   if ( $debug ) echo "++ creating spec file '$spec_file'..."

   (echo "# delimits comments" > $spec_file) >& /dev/null

   if ( $status ) then
       set endstr = ( "failure: no permissions to create spec file" \
          "'$suma_dir/$spec_file'" )
       goto L_BAD_END
   endif

   # note user, date, machine, pwd, command line
   echo ""                     >> $spec_file
   echo "# Creation information:"            >> $spec_file
   echo "#     user    : $user"            >> $spec_file
   echo "#     date    : `date`"            >> $spec_file
   echo "#     machine : `uname -n`"         >> $spec_file
   echo "#     pwd     : $cwd"            >> $spec_file
   echo "#     command : $prog_name $argv"         >> $spec_file
   echo ""                     >> $spec_file
   
   # define the group
   echo "# define the group"            >> $spec_file
   echo "        Group = $subj_id"            >> $spec_file
   echo ""                     >> $spec_file

   # define the states
   echo "# define various States"            >> $spec_file
   foreach attrib ( $surf_attribs )
      if (  $attrib =~ *sphere*  || \
            $attrib =~ *flat*    || \
            $attrib =~ *infla*) then
         echo "        StateDef = ${attrib}_lh"         >> $spec_file
         echo "        StateDef = ${attrib}_rh"         >> $spec_file
      else
         echo "        StateDef = $attrib"         >> $spec_file
      endif
   end
   echo ""                     >> $spec_file
   foreach hand ( lh rh )
   if ( $hand == lh ) then
       set list_cur = ( $list_lh )      # get a current list copy
   else
       set list_cur = ( $list_rh )
   endif

   foreach surf ( $list_cur )
       set s_head = `echo $surf | cut -d. -f1`   # up to first '.'
       set s_state = `echo $surf | cut -d. -f2-`   # after first '.'
       
       if ( $s_state =~ *sphere* || \
            $s_state =~ *flat*   || \
            $s_state =~ *infla*) then
         set s_state = "${s_state}_${hand}"
       endif
       
      set label_dset = ""
       # check for SAME mapping ref
       if ( "$surf.${sex}" == "$s_head.smoothwm.${sex}" ) then
      set map_ref = SAME
         if ( -f ${hand}.aparc.a2005s.annot.niml.dset ) then
      set label_dset = ${hand}.aparc.a2005s.annot.niml.dset
         endif
         if ( -f ${hand}.aparc.a2009s.annot.niml.dset ) then
      set label_dset = ${hand}.aparc.a2009s.annot.niml.dset
         endif
       else
      set map_ref = $s_head.smoothwm.${sex}
       endif

       # EmbedDimension is 2 for .flat surfaces, else 3
       if ( $surf =~ *.flat* ) then
      set embed_ref = 2
       else
      set embed_ref = 3
       endif

       # Anatomical flag
       if ( $surf =~ *.white* ||\
            $surf =~ *.smoothwm* ||\
            $surf =~ *.pial*) then
      set anatomical = Y
       else
      set anatomical = N
       endif

       echo "NewSurface"               >> $spec_file
       echo "        SurfaceFormat = ASCII"      >> $spec_file
       echo "        SurfaceType = ${stp}"      >> $spec_file
       echo "        ${sfieldname} = $surf.${sex}"   >> $spec_file
       echo "        LocalDomainParent = $map_ref"      >> $spec_file
       echo "        SurfaceState = $s_state"      >> $spec_file
       echo "        EmbedDimension = $embed_ref"      >> $spec_file
       echo "        Anatomical = $anatomical" >> $spec_file
         if ("$label_dset" != "") \
       echo "        LabelDset = $label_dset" >> $spec_file
       echo ""                  >> $spec_file
   end
   end    # foreach hand

   echo "++ created spec file'$suma_dir/$spec_file'"


    cd $start_dir

goto L_CREATE_SPEC_DONE


####################################################################
# actually create an AFNI dataset

L_CREATE_BRICK:

    cd $orig_dir
    # verify existence of 256 COR files, Checking for 256 is now obsolete ZSS Oct 14 04
    set cor_files = ( `echo COR-???` ) >& /dev/null 
    
    if ( $status || $use_mgz == 1) then
      if ($orig_mgz != "") then   
         #Try the orig.mgz business
         echo "No COR files, converting $orig_mgz ..."
         cd $start_dir
         mri_convert -ot cor $orig_mgz $orig_dir
         cd $orig_dir
      else 
         if ($use_mgz == 0) then
            set endstr = "failure: No COR files and no $orig.mgz"
         else
            set endstr = "failure: No mgz files"
         endif
         goto L_BAD_END
      endif
    else
      if ($orig_mgz != "" && $neuro_ori != -1) then  
         echo ""
         echo "Notice: Found COR files and $orig_mgz"
         echo "******** Will use COR files."
         echo ""
      else if ($orig_mgz == "" && $neuro_ori == -1) then  
         set endstr = "No orig_mgz found for -nocor option"
         goto L_BAD_END
      endif
    endif
   
   #check again
   set cor_files = ( `echo COR-???` ) >& /dev/null 
    if ( $status ) then
   set endstr = "failure: did not find COR files under '$orig_dir'"
   goto L_BAD_END
    #else if ( $#cor_files != 256 ) then
   #set endstr = "failure: 256 COR files required under '$orig_dir'"
   #goto L_BAD_END
    endif

    # create BRIK with to3d

    #SLAB changed by ZSS Mar 12 04. was 0.5 mm off. Bug pointed
    #out by A. Thomas. Good info found here:
    #www.wideman-one.com/gw/brain/fs/coords/fscoords.htm
    
      #New block added for reading COR images that are no longer 
      #256x256x256 / 1mm cubic voxels ZSS Oct 20 04
      
      set Sr = `grep -w imnr0 COR-.info`
      set i0 = $Sr[$#Sr]
      set Sr = `grep -w imnr1 COR-.info`
      set i1 = $Sr[$#Sr]
      set Sr = `grep -w fov COR-.info`
      set FOV = `ccalc -eval $Sr[$#Sr] \* 1000`
      set Sr = `grep -w x COR-.info | grep -v xform`
      set nx = $Sr[$#Sr]
      set Sr = `grep -w y COR-.info`
      set ny = $Sr[$#Sr]
      set Sr = `grep -w thick COR-.info`
      set zthick = `ccalc -eval $Sr[$#Sr] \* 1000`
      set Sr = `grep -w psiz COR-.info`
      set psize = `ccalc -eval $Sr[$#Sr] \* 1000`

      set xSLAB0 = `ccalc -eval $nx \* $psize / 2`
      set xSLAB1 = `ccalc -eval $xSLAB0 - $psize`
      set ySLAB0 = `ccalc -eval $ny \* $psize / 2`
      set ySLAB1 = `ccalc -eval $ySLAB0 - $psize`
      set zSLAB0 = `ccalc -eval '('$i1 - $i0 + 1')' \* $zthick / 2`
      set zSLAB1 = `ccalc -eval $zSLAB0 - $zthick`

    
    if ( $neuro_ori == 1 ) then
      to3d -prefix $afni_prefix \
               -xSLAB $xSLAB0'L'-$xSLAB1'R' \
               -ySLAB $ySLAB0'S'-$ySLAB1'I' \
               -zSLAB $zSLAB0'P'-$zSLAB1'A' '3Db:0:0:'$nx':'$ny':1:COR-???'
      #to3d -prefix $afni_prefix -xSLAB 127L-128R -ySLAB 128S-127I \
      #  -zSLAB 128P-127A $cor_files
    else if ( $neuro_ori == 0 ) then
      to3d -prefix $afni_prefix \
               -xSLAB $xSLAB0'R'-$xSLAB1'L' \
               -ySLAB $ySLAB0'S'-$ySLAB1'I' \
               -zSLAB $zSLAB0'P'-$zSLAB1'A' '3Db:0:0:'$nx':'$ny':1:COR-???'
      #to3d -prefix $afni_prefix -xSLAB 128R-127L -ySLAB 128S-127I \
      #  -zSLAB 128P-127A $cor_files
    else if ( $neuro_ori == -1 ) then
      if ( "$orig_mgz" == '' ) then
         set endstr = "Could not locate orig_mgz volume"
         goto L_BAD_END
      endif
      mri_convert -ot nii $start_dir/$orig_mgz \
                  $start_dir/$suma_dir/${orig_mgz:t:r}".nii"
    else 
      set endstr = "Bad value for neuro_ori of $neuro_ori"
      goto L_BAD_END
    endif

    if ( $neuro_ori == -1 ) then
       mv $start_dir/$suma_dir/${orig_mgz:t:r}".nii"  \
            $start_dir/$suma_dir/$afni_dataset.nii
       echo "++ created AFNI file '$suma_dir/$afni_dataset.nii'"
    else
       mv $afni_dataset.HEAD $afni_dataset.BRIK* $start_dir/$suma_dir

       echo "++ created AFNI file '$suma_dir/$afni_dataset.HEAD'"
       echo "++ created AFNI file '$suma_dir/$afni_dataset.BRIK'"
    endif

    cd $start_dir
    
    if ($#other_mgz > 0) then
       #Do the other .mgz files
       foreach other ($other_mgz)
         #@SUMA_FSvolToBRIK $other $start_dir/$suma_dir/$other:t:r
         #rm -f $start_dir/$suma_dir/COR-??? $start_dir/$suma_dir/COR-.info \
         #                                                        >& /dev/null
         mri_convert -ot nii $other $start_dir/$suma_dir/${other:t:r}".nii"
         #Now shift the damned volumes to match the surface volumes
         cd $start_dir/$suma_dir/
         if ($neuro_ori != -1 ) then
            @Align_Centers -base $afni_dataset -dset ${other:t:r}".nii" -no_cp
         endif
         if ( "${other:t:r}" == "aparc.a2005s+aseg" ||   \
              "${other:t:r}" == "aparc.a2009s+aseg" ||   \
              "${other:t:r}" == "aseg" ||   \
              "${other:t:r}" == "aparc+aseg" ) then
            3dmerge -1rank -prefix "${other:t:r}_rank.nii" ${other:t:r}".nii"
            if (0) then
               # kill this block soon.
               set kk = `echo ${other:t:r} | cut -f 1 -d '+'`
               mv ${kk}.rankmap.1D ${other:t:r}.rankmap.1D 
            endif
            @FS_roi_label  -name ALL \
                           -rankmap ${other:t:r}_rank.rankmap.1D \
                           -labeltable ${other:t:r}_rank \
                              > ${other:t:r}_rank.niml.lt.log
            if ($status) then
               echo "Failed to make roi labels for $other"
            else
               3drefit -labeltable ${other:t:r}_rank.niml.lt \
                                  "${other:t:r}_rank.nii"
               @MakeLabelTable -atlasize_labeled_dset "${other:t:r}_rank.nii"
            endif
         endif
         cd - 
       end
    endif
    
    #Change from orig space
    if ( ${set_space} != '') then
      cd $start_dir/$suma_dir/
      3drefit -space ${set_space} *.nii* *.HEAD
      cd -
    endif
      
goto L_CREATE_BRICK_DONE

####################################################################
# echo run MapIcosahedron 

L_ICO:

   if ($#ldlist < 1) goto L_ICO_DONE
   
   echo ""
   echo "------------------------------------------------------------------"
   echo "Running MapIcosahedron for both hemispheres at ld values of $ldlist"
   cd $start_dir/$suma_dir/
   set ild = 0
   foreach ld ($ldlist)
      @ ild ++
      set both = ()
      foreach hem (lh rh)
         set sp = {$subj_id}_${hem}.spec
         set dm = ()
         if ($#fsdsets > 0) then
            foreach fsd ($fsdsets)
               if ( -f ${hem}.${fsd}.gii.dset ) then
                  set dm = ($dm "-dset_map" ${hem}.${fsd}.gii.dset)
               endif
            end
         endif
         if ( -f $sp) then
            MapIcosahedron -spec $sp -ld $ld $dm -prefix $ldpref[$ild] 
            set both = ($both ${ldpref[$ild]}{$subj_id}_${hem}.spec)
         endif
      end
      if ($#both == 2) then
         inspec -LRmerge $both -prefix ${ldpref[$ild]}{$subj_id}_both.spec
      endif
   end
   cd -
   
goto L_ICO_DONE

####################################################################
# echo details for the user to launch suma and afni, in order to
# check the alignment

L_TEST_SURF_VOL:

    echo ""
    echo "------------------------------------------------------------------"
    echo "Please verify that the datasets are aligned properly in both"
    echo "afni and suma.  You may do this by running the following commands:"
    echo ""
    echo "    cd $suma_dir"
    echo "    afni -niml &"
    if ( -f ${suma_dir}/${afni_dataset}.nii ) then
      set ovol = ${afni_dataset}.nii
    else if ( -f ${suma_dir}/${afni_dataset}.nii.gz ) then
      set ovol = ${afni_dataset}.nii.gz
    else 
      set ovol = $afni_dataset
    endif
    echo "    suma -spec $spec_file -sv $ovol"

goto L_TEST_SURF_VOL_DONE


####################################################################
# display help and exit

L_HELP_END:

echo ""
echo "$prog_name - prepare for surface viewing in SUMA"
echo ""
echo "    This script goes through the following steps:"
echo "      - verify existence of necessary programs "
echo "        (afni, to3d, suma, mris_convert)"
echo "      - determine the location of surface and COR files"
echo "      - creation of ascii surface files via 'mris_convert'"
echo "      - creation of left and right hemisphere SUMA spec files"
echo "      - creation of an AFNI dataset from the COR files via 'to3d'"
echo "      - creation of AFNI datasets from various .mgz volumes created"
echo "        by FreeSurfer. The segmentation volumes with aseg in the "
echo "        name are best viewed in AFNI with the FreeSurfer_Seg_255"
echo "        colormap. See bottom of @SUMA_FSvolToBRIK -help for more"
echo "        info." 
echo ""
echo "      - all created files are stored in a new SUMA directory"
echo ""
echo "  Usage: $prog_name [options] -sid SUBJECT_ID"
echo ""
echo "  examples:"
echo ""
echo "    $prog_name -sid subject1"
echo "    $prog_name -help"
echo "    $prog_name -fspath subject1/surface_stuff -sid subject1"
echo "    $prog_name -neuro -sid 3.14159265 -debug 1"
echo ""
echo "  options:"
echo ""
echo "    -help    : show this help information"
echo ""
echo "    -debug LEVEL    : print debug information along the way"
echo "          e.g. -debug 1"
echo "          the default level is 0, max is 2"
echo ""
echo "    -fspath PATH    : path to 'surf' and 'orig' directories"
echo "          e.g. -fspath subject1/surface_info"
echo "          the default PATH value is './', the current directory"
echo ""
echo "          This is generally the location of the 'surf' directory,"
echo "          though having PATH end in surf is OK.  The mri/orig"
echo "          directory should also be located here."
echo ""
echo "          Note: when this option is provided, all file/path"
echo "          messages will be with respect to this directory."
echo ""
echo "    -use_mgz        : use MGZ volumes even if COR volumes are there"
echo ""
echo "    -neuro          : use neurological orientation"
echo "          e.g. -neuro"
echo "          the default is radiological orientation"
echo ""
echo "          In the default radiological orientation, the subject's"
echo "          right is on the left side of the image.  In the"
echo "          neurological orientation, left is really left."
echo ""
echo "    -nocor: This option is no longer supported because it created"
echo "            GIFTI surfaces with coordinates in RAI, rather than LPI"
echo "            which is the GIFTI standard. While using RAI surfaces"
echo "            within AFNI/SUMA is not problematic, the resultant GIFTI"
echo "            surfaces do not port well to other software."
echo "            The replacement option for -nocor is -GNIFTI but the"
echo "            surfaces will have negated coordinates along the x and y"
echo "            compared to those with -nocor."
echo "            GIFTI surfaces produced with SUMA programs compiled before"
echo "            August 1st 2013 will have their X and Y coordinates "
echo "            negated and will no longer line up with the anatomy. "
echo "            Correcting such surfaces can be done with ConvertSurface"
echo "            with the following command:"
echo "              ConvertSurface  -i lh.smoothwm.gii -o_gii lh.smoothwm \"
echo "                              -overwrite -xmat_1D NegXY"
echo "            or for an entire SUMA directory:"
echo "              cd SUMA"
echo "              tcsh"
echo '              foreach ss (*.gii)'
echo '                 ConvertSurface  -i $ss -o_gii $ss \'
echo "                                 -overwrite -xmat_1D NegXY"
echo "              end"
echo ""            
echo "    -GNIFTI:Produce files in exchangeable formats. With this option"
echo "    -NIFTI :COR volumes are no longer used and output volumes"
echo "    -GIFTI :and surfaces are in alignment with the original "
echo "    -IFTI  :volume used to create the surface. All volumes are"
echo "            written out NIFTI format, and all surfaces are"
echo "            in GIFTI format."
echo "            This option is incompatible with -neuro or -use_mgz"
echo "    NOTE for -GNIFTI:"
echo "        If you really care that the volumes in SUMA/ are in exact"
echo "        register with the volume you passed to FreeSurfer, you "
echo "        should be sure that the volume passed to FreeSurfer has "
echo "        an even number of slices in all directions and a voxel"
echo "        resolution of 1x1x1, otherwise the resultant volumes in"
echo "        SUMA/ might be off by half a voxel or less in directions"
echo "        with odd number of slices. The reason has to do (I think)"
echo "        with FreeSurfer's resampling and volume centering approach."
echo "        In either case, surfaces and volumes under SUMA/ will be in"
echo "        proper register."
echo "        For example, when creating a surface model of the TT_N27 brain"
echo "        I padded the TT_N27+tlrc volume before submitting it to "
echo "        FreeSurfer with the following command:"
echo "           3dZeropad -L 1 -P 1 -S 1 -prefix anat.nii TT_N27+tlrc.HEAD"
echo "        After zeropadding, anat.nii remains in perfect register with"
echo "        TT_N27+tlrc by it has an even number of slices in all"
echo "        directions: 3dinfo -n4  -d3 -prefix anat.nii"
echo "          162  192  152   1    1.0     1.0     1.0     anat.nii"
echo ""
echo "    -inflate INF: Create modereately inflated surfaces using"
echo "                  SurfSmooth. INF controls the amount of smoothness"
echo "                  in the final image. It is the number of iterations"
echo "                  in the command such as: "
echo "             SurfSmooth  -i lh.white.asc    -met NN_geom \"
echo "                    -Niter 200  -o_gii  -surf_out lh.inf_200 \"
echo "                    -match_vol 0.01"
echo "                  You can use multiple instances of -inflate to create"
echo "                  inflations of various levels."
echo "    -set_space  SPACE: Set the space flag of all volumes to"
echo "            SPACE (orig, MNI, TLRC, MNIa). The default is "
echo "            orig space."
echo "            You should only use this option when the volume you"
echo "            passed to FreeSurfer was not in 'orig' space."
echo "            Use '3dinfo -space YOUR_DATASET' to find the space "
echo "            of a certain dataset." 
echo ""
echo "    -sid SUBJECT_ID : required subject ID for file naming"
echo ""
echo "    -ld LD : Create standard mesh surfaces with mesh density"
echo "             linear depth (see MapIcosahedron -help, option -ld)"
echo "             set to LD. You can use multiple -ld options."
echo "             By default the script will run ld values of 141 and"
echo "             60."
echo "    -ldpref LDpref: Supply what ends up being the -prefix option"
echo "                    for MapIcosahedron. By default it is std.LD."
echo "                    You need as many -ldpref as you have -ld"
echo "    -no_ld: Do not run MapIcosahedron."
echo ""
echo "  Making use of FreeSurfer's -contrasurfreg output with MapIcosahedron:"
echo "  This script will create SUMA versions of lh.rh.sphere.reg and "
echo "  rh.lh.sphere.reg but in this current state, MapIcosahedron does"
echo "  not attempt to use them for backward compatibility."
echo "  Should you want to create standard mesh surfaces with node"
echo "  index correspondence across the hemispheres you will need to run"
echo "  MapIcosahedron manually in the output SUMA/ directory. \n"
echo "  For example:"
echo "      MapIcosahedron    -spec SUBJ_rh.spec -ld 60 \"
echo "                        -dset_map rh.thickness.gii.dset \"
echo "                        -dset_map rh.curv.gii.dset \"
echo "                        -dset_map rh.sulc.gii.dset \"
echo "                        -morph rh.lh.sphere.reg.gii \"
echo "                        -prefix std.60.lhreg."
echo "      This command is very similar to the one use to create the"
echo "      default output spec file std.60.SUBJ_rh.spec (look at the "
echo "      top of the spec file for a record of the command that created it),"
echo "      except for the last two options -morph and -prefix."
echo "      By using -morph rh.lh.sphere.reg.gii the resultant standard-mesh "
echo "      right hemispheres (std.60.lhreg.rh.*.gii) will have node index "
echo "      correspondence with std.60.lh.*.gii surfaces."
echo "      To verify visually the correspondence, run the following:"
echo "         count -column 0 36001 > std.60.lh.rh.nodeindex.1D.dset"
echo '         suma -noniml -spec std.60.SUBJ_lh.spec &'
echo '         suma -noniml -spec std.60.SUBJ_rh.spec &'
echo '         suma -noniml -spec std.60.lhreg.SUBJ_rh.spec &'
echo '       Then load std.60.lh.rh.nodeindex.1D.dset into each of the three'
echo "       SUMA windows. Note how the color pattern (node indices) matches"
echo "       between SUBJ_lh and lhreg.SUBJ_rh surfaces, "
echo "       but NOT between SUBJ_lh and SUBJ_rh surfaces."
echo ""
echo "  notes:"
echo ""
echo "    0. More help may be found at https://afni.nimh.nih.gov/ssc/ziad/SUMA/SUMA_doc.htm"
echo "    1. Surface file names should look like 'lh.smoothwm'."
echo "    2. Patches of surfaces need the word patch in their name, in"
echo "       order to use the correct option for 'mris_convert'."
echo "    3. Flat surfaces must have .flat in their name."
echo "    4. You can tailor the script to your needs. Just make sure you rename it or risk"
echo "       having your modifications overwritten with the next SUMA version you install."
echo ""
echo "     R. Reynolds (rickr@codon.nih.gov)"
echo "     Z. Saad (saadz@mail.nih.gov)"
echo "     M. Beauchamp (Michael.S.Beauchamp@uth.tmc.edu)"
echo ""

    exit


####################################################################
# failure!

L_BAD_END:

    echo ""
    if ( "$endstr" != "" ) echo "$endstr"
    echo program failure: exiting...
    echo ""

    exit

####################################################################
# finished!

L_GOOD_END:

    if ( $debug > 1 ) unset echo

    echo ""
    if ( "$endstr" != "" ) then
   echo "$endstr"
   echo ""
    endif

    exit   # just to be consistent