#!/bin/tcsh -f

PARSE: 
   echo "Parsing ..."
   set Narg = $#
   
   set tmp = ___tmp_WO_
   set pd_default = 0
   set cnt = 1
   set ref_in = ''
   set anat_in = ''
   set KeepTmp = 0
   set verb = 0
   set pd = $pd_default
   set suf = ""
   set ss_opt = ''
   set dxyz_list = ( 1 1 1 )
   
   if ("$1" == '') goto HELP

   while ($cnt <= $Narg)
      set donext = 1;
      
      if ($donext && "$argv[$cnt]" == "-help" || "$argv[$cnt]" == "-h") then
         goto HELP
      endif
      
      if ($donext && "$argv[$cnt]" == "-dxyz") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need voxel size after -dxyz"
            goto END
         else
            @ cnt ++
            set dxyz_list[1] = "$argv[$cnt]"
            set dxyz_list[2] = "$argv[$cnt]"
            set dxyz_list[3] = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      
      if ($donext && "$argv[$cnt]" == "-dx") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need voxel size after -dx"
            goto END
         else
            @ cnt ++
            set dxyz_list[1] = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      if ($donext && "$argv[$cnt]" == "-dy") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need voxel size after -dy"
            goto END
         else
            @ cnt ++
            set dxyz_list[2] = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      if ($donext && "$argv[$cnt]" == "-dz") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need voxel size after -dz"
            goto END
         else
            @ cnt ++
            set dxyz_list[3] = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      
      if ($donext && "$argv[$cnt]" == "-suffix") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need a string after -suffix"
            goto END
         else
            @ cnt ++
            set suf = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      
      if ($donext && "$argv[$cnt]" == "-base") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need template volume after -base"
            goto END
         else
            @ cnt ++
            set ref_in = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif

      if ($donext && ( "$argv[$cnt]" == "-pad_base" || "$argv[$cnt]" == "-pad_input" ) ) then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need value after -pad_base (or -pad_input) "
            goto END
         else
            @ cnt ++
            set pd = "$argv[$cnt]"
            if ($pd > 100 || $pd < 0) then
               echo "Error: -base_pad should be between 0 and 100"
               echo "I have $pd"
               goto END
            endif
            set donext = 0   
         endif   
      endif
      if ($donext && "$argv[$cnt]" == "-input") then
         set pLoc = $cnt      
         if ($pLoc == $Narg) then
            echo "Need  volume after -input"
            goto END
         else
            @ cnt ++
            set anat_in = "$argv[$cnt]"
            set donext = 0   
         endif   
      endif
      if ($donext && "$argv[$cnt]" == "-keep_tmp") then
         set KeepTmp = 1;      
         set donext = 0      
      endif
      if ($donext && "$argv[$cnt]" == "-clean") then
         rm -f ${tmp}_*.* >& /dev/null
         goto END      
         set donext = 0      
      endif
      if ($donext && "$argv[$cnt]" == "-verb") then
         set verb = 1;      
         set donext = 0      
      endif
      if ($donext == 1) then
         echo "Error: Option or parameter '$argv[$cnt]' not understood"
         goto END
      endif
      @ cnt ++
   end
   if ("$anat_in" == "" || "$ref_in" == "") then
      echo "Error: Need both -base and -input parameters"
      goto END
   endif

set base_parsed = `@parse_afni_name $ref_in`
set anat_parsed = `@parse_afni_name $anat_in`

set base_o = $base_parsed[2]
set anat_o = $anat_parsed[2]

set base = $base_o
set anat = $anat_o

goto MOVE_BASE

ALIGN_CENTERS:
   #Not needed, assuming there are no clipping problems.
   #Better do this for all oblique sets outside of this script if necessary
   rm -f ${anat}_shft*.???? >& /dev/null
   @Align_Centers -base ${base}+orig -dset ${anat}+orig
   set anat = ${anat}_shft

MOVE_BASE:
#Try different large angles, moving base because it is much larger
   rm -f ${tmp}_${base}_rs*.???? >& /dev/null
   3dresample -dxyz $dxyz_list[1] $dxyz_list[2] $dxyz_list[3]   -inset ${base_parsed[1]}/${base}+orig -prefix ./${tmp}_${base}_rs
   set base = ${tmp}_${base}_rs
   rm -f ${tmp}_${anat}_rs*.???? >& /dev/null
   3dresample -inset ${anat_parsed[1]}/${anat}+orig -prefix ./${tmp}_${anat}_rs -master ${base}+orig
   set anat = ${tmp}_${anat}_rs
   set anat_i = $anat
   set base_i = $base
   set olst = ''
   set ilst = ''
   set RMSmin = 9000;
   set ANGmin = -1;
   set CNTmin = 0;
   set cnt = 0
   foreach ang ( 0 20 40 )
      echo "Trying ${ang}R 0A 0I...(`date`)"
      @ cnt ++
      set anat = ${anat_i}
      rm -f ${base}_R${ang}*.???? >& /dev/null
      3drotate -rotate ${ang}R 0A 0I -prefix ${base}_R${ang} ${base}+orig
      set ilst = ($ilst ${base}_R${ang}+orig)
      rm -f ${anat}_b${ang}_aa*.???? >& /dev/null
         
         3dvolreg -twopass -cubic -prefix ${anat}_b${ang}_aa -dfile ${anat}_b${ang}_aa.log -base ${base}_R${ang}+orig ${anat}+orig.
         #get the rms result 
         set dfile = `cat ${anat}_b${ang}_aa.log`
         if ( `ccalc -form int -expr step"($RMSmin - $dfile[9])"` ) then
            set ANGmin = $ang;
            set RMSmin = $dfile[9];
            set CNTmin = $cnt;
         endif
         set anat = ${anat}_b${ang}_aa
      
      set olst = ($olst ${anat}+orig)
   end
   echo $olst
   echo $ilst
   echo "`date`:Optimal registration found at test number $CNTmin ANGmin = $ANGmin, RMSmin = $RMSmin" >> ${tmp}_LeLog.txt
   tail ${tmp}_LeLog.txt

goto THE_FINAL_GIG

STORE_RESULTS:
   #Not necessary unless you are debugging
   rm -f ${anat_i}_aa_all*.???? ${base_i}_all*.???? >& /dev/null
   3dTcat -prefix ${anat_i}_aa_all $olst
   3dTcat -prefix ${base_i}_all $ilst

   #afni -noplugouts ${base}+orig.HEAD ${anat_i}_aa_all+orig.HEAD ${base_i}_all+orig.HEAD   &

APPLY_BEST_MOVE_BASE:
   #get inverse transform of best alignment
   cat_matvec "$olst[$CNTmin]"::VOLREG_MATVEC_000000 -I > ${tmp}_iVR_MV.1D
   #apply it to the pre-moved base
   3drotate -matvec_dicom ${tmp}_iVR_MV.1D -prefix `@GetAfniPrefix $ilst[$CNTmin]`_alnd_${anat_i} $ilst[$CNTmin]

   echo "Run: afni -noplugouts `@GetAfniPrefix $ilst[$CNTmin]`_alnd_${anat_i}+orig.HEAD ${anat_i}+orig.HEAD"

THE_FINAL_GIG:
   if ("$suf" == "") then
      set suf = "_alnd_${anat_o}"
   endif

   #apply full xform to initial base
   cat_matvec ${base}_R${ANGmin}+orig::ROTATE_MATVEC_000000 "$olst[$CNTmin]"::VOLREG_MATVEC_000000 -I > ${tmp}_BaseToAnat.1D
   #apply to initial base
   rm -f ${base_o}${suf}*.???? >& /dev/null
   3drotate -matvec_dicom ${tmp}_BaseToAnat.1D -prefix ./${base_o}${suf} ${base_parsed[1]}/${base_o}+orig 
   
   echo "Run: afni -noplugouts ${base_o}_alnd_${anat_o}+orig.HEAD ${anat_o}+orig.HEAD"

goto CLEANUP


FINE_REG:
#Cannot do finer registration because scaling is not allowed, partial coverage...
if (0) then
      3dWarpDrive -coarserot -base `@GetAfniPrefix $ilst[$CNTmin]`_alnd_${anat_i}+orig\
               -prefix Anat_alnd_MovedBase -affine_general \
               -input ${anat_i}+orig -quintic -verb
endif

CLEANUP:
if ($KeepTmp == 0) then
   rm -f ${tmp}_*.* >& /dev/null
endif

goto END

HELP:
   if (1) then
      echo "Usage 1: A script to align a full coverage T1 weighted non-oblique dataset"
      echo "         to match a partial coverage T1 weighted non-oblique dataset "
      echo "         Alignment is done with a rotation and shift (6 parameters) transform"
      echo "         only."
      echo ""
      echo " Script is still in testing phase"
      echo "" 
      echo "   `basename $0` [options] <-base FullCoverageT1> <-input PartialCoverageObliqueT1>"
      echo "   Mandatory parameters:"
      echo "      -base  FullCoverageT1:  Reference anatomical full coverage volume."
      echo ""
      echo "      -input  PartialCoverageObliqueT1:  The name says it all."
      echo ""                        
      echo "   Optional parameters:"
      #echo "      -pad_base  MM  :  Pad the base dset by MM mm in each directions."
      #echo "                        That is needed to  make sure that datasets"
      #echo "                        requiring wild rotations do not get cropped."
      #echo "                        Default is MM = $pd_default."
      #echo "                        If your output dataset is clipped, try increasing"
      #echo "                        MM to `ccalc -eval $pd_default+10` or `ccalc -eval $pd_default+20`."
      echo "      -suffix  SUF   :  Output dataset name is formed by adding SUF to"
      echo "                        the prefix of the base dataset."
      echo "                        The default suffix is _alnd_PartialCoverageObliqueT1"
      echo "      -keep_tmp      :  Keep temporary files."
      echo "      -clean         :  Clean all temp files, likely left from -keep_tmp"
      echo "                        option then exit."
      echo "      -dxyz MM          : Cubic voxel size of output DSET in TLRC"
      echo "                          space Default MM is 1. If you do not"
      echo "                          want your output voxels to be cubic"
      echo "                          Then use the -dx, -dy, -dz options below."
      echo "      -dx MX            : Size of voxel in the x direction"
      echo "                          (Right-Left). Default is 1mm."
      echo "      -dy MY            : Size of voxel in the y direction"
      echo "                          (Anterior-Posterior). Default is 1mm."
      echo "      -dz MZ            : Size of voxel in the z direction."
      echo "                          (Inferior-Superior).Default is 1mm."      
      echo "   Example:"
      echo "   `basename $0` -base ah_SurfVol+orig. -input ah_T1W_anat+orig."
      echo ""
      echo ""
      echo "Written by Ziad S. Saad, for Ikuko (ziad@nih.gov)"
      echo "                        SSCC/NIMH/NIH/DHHS"
      echo ""
      goto END
   endif


END: