#!/bin/csh -f
#
# ggfm
#
# Pre-process GE DICOM files and then calls epidwarp4.ucsd to perform unwarping
# 
#
# Requires: FSL tools from fMRIB 
#           AFNI tools from NIH
#           dicomrx from cfmriweb.ucsd.edu/fmap/dicomrx
#
# Calls:  epidewarp4.ucsd and dicomrx
#
# What this script does:
#  1) Creates NIFTI from DICOM files. 
#  2) Determines echo times from DICOM headers
#  3) Converts field map DICOM to NIFTI
#  4) Generates complex AVW volumes and magnitudes for masking
#  5) Determine EPI dwell time from DICOM header
#  6) Converts EPI DICOM files to NIFTI
#  7) Calls epidewarp4.ucsd to unwarp the files
#
# New capabilities:
#  1. multiple coils
#
# Version History
#  1.0 04/19/2010 KL created by modifying the ppge4 code to work with the fm_grass protocol. 
#  1.1 06/17/2010 KL added -n option to allow users to supply their own epivol in .nii.gz format.
# Send Comments/Questions to kunlu@ucsd.edu.
#

set VERSION = '$Id$';
set inputargs = ($argv);
set DateStr = "`date '+%y%m%d%H%M'`"

set fm = ();
#set dir1 = (); 
#set dir2 = (); 
set instem = ();
set tediff = ();
set tmpdir = 'tmp';
set outstem = ();
set TR = ();
set merge_in = 1;
set domask = 1;
set dobrik = 0;
set unwarpdir = ();
set domoco = 1;
set refnum = ();
set docleanup = 1;
set epinifti = (); 
set usrnifti = 0;

set ncl=-1;
set nphase = 8;
set betf=0.5;
set reps = 0;
set nslices = 0;
set cleanup = 1;
set PrintHelp = 0;
set dopmask = 1;
set dofmask = 1;
set dounwrap = 1;
set dogzip = 1;
set fmap = ();
set postfix = '.nii.gz';

if($#argv == 0) goto usage_exit;
set n = `echo $argv | grep version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif
set n = `echo $argv | grep help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
endif

goto parse_args;
parse_args_return:


goto check_params;
check_params_return:


## Get/Create tmp directory ##
mkdir -p $tmpdir
if( ! -e $tmpdir) then
  echo "ERROR: cannot find tmp dir $tmpdir"
  exit 1;
endif 
# rm -rf $tmpdir/*


set curdir = `pwd`;

echo " "
echo " Grabbing TE1 and TE2 Data"
echo ""

# convert to NIFTI format
echo $fm

# determine # field map files
set nfiles = `ls -1 $curdir/$fm/i* | wc -l`
echo "$nfiles images in total"


#determine the type of recerivers used
#Changed from dicom2 to dicom_hdr (afni) for the rest of us.
# -cbk
    set f1 = `ls $fm/i*.1`
    set l1 = `dicom_hdr $f1 | grep "0018 1250"`
    set coil  = `echo "$l1" | awk -F// '{print $3}'`
   
    echo "INFO: $coil receivers was used to acquire the field map" 
   
@ cntd = 6

# assuming mag,ph,re,im
set re_te1 = re_te1$postfix
set im_te1 = im_te1$postfix
set re_te2 = re_te2$postfix
set im_te2 = im_te2$postfix

@ cnt0 = 2
cd $curdir/$fm/
set cmd= (to3d -prefix $re_te1 `count -digits 1 -root 'i*.' $cnt0 $nfiles $cntd`)
# echo $cmd 
$cmd
if($status) exit 1;
@ cnt0 = 3
set cmd= (to3d -prefix $im_te1 `count -digits 1 -root 'i*.' $cnt0 $nfiles $cntd`)
# echo $cmd 
$cmd
if($status) exit 1;
@ cnt0 = 5
set cmd=(to3d -prefix $re_te2 `count -digits 1 -root 'i*.' $cnt0 $nfiles $cntd`)
# echo $cmd 
$cmd 
if($status) exit 1;
@ cnt0 = 6
set cmd=(to3d -prefix $im_te2 `count -digits 1 -root 'i*.' $cnt0 $nfiles $cntd`)
# echo $cmd 
$cmd 
if($status) exit 1;

mv *$postfix $curdir/$tmpdir


# make magnitude images of te1 for masking purpose later.
cd $curdir/$fm/
set mag_te1 = mag_te1$postfix

@ cnt0 =  + 1

set cmd=(to3d -prefix $mag_te1 `count -digits 1 -root 'i*.' $cnt0 $nfiles $cntd`)
echo $cmd
$cmd 
if($status) exit 1;
mv $mag_te1 $curdir/$tmpdir

cd $curdir


# determine echo times
# Changed dicom2 to dicom_hdr
# get tediff from opuser2 KL
if ( $tediff == "auto" ) then
    set f1 = `ls $fm/i*.1`
    set l1 = `dicom_hdr $f1 | grep "0019 10a8"`
    set te1  = `echo "$l1" | awk -F// '{print $3}'`
    set tediff = ` echo " 1000.0*($te1) " | bc -l`; 
endif
echo "INFO: tediff = $tediff usec" 

# calculate conjugate product

echo "calculate conjugate product"
3dcalc -a $tmpdir/re_te2$postfix -b $tmpdir/im_te2$postfix -c $tmpdir/re_te1$postfix -d $tmpdir/im_te1$postfix -expr "(a*c+b*d)/1000" -prefix $tmpdir/re_dph$postfix


3dcalc -a $tmpdir/re_te2$postfix -b $tmpdir/im_te2$postfix -c $tmpdir/re_te1$postfix -d $tmpdir/im_te1$postfix -expr "(-a*d+b*c)/1000" -prefix $tmpdir/im_dph$postfix


# generate complex AVW volumes of the conjugate product
set cp_dph = $tmpdir/cp_dph$postfix
fslcomplex -complex  $tmpdir/re_dph$postfix $tmpdir/im_dph$postfix $cp_dph


#Also make up an AFNI brik to act as a geometric parent for to3d later on
set nslices = `fslinfo $tmpdir/$mag_te1 | awk '{if($1 == "dim3") print $2}'`
set gp = $tmpdir/$mag_te1

#determine EPI dwell time
#Changed dicom2 to dicom_hdr
set f1 = `ls $epidir/*.1 `
set l1 = `dicom_hdr $f1 | grep -i "0043 102c"`
set dwell  = `echo "$l1" | awk -F// '{print $3}'`
echo ""
echo "INFO: EPI dwell time = $dwell usec"

set asset = `dicom_hdr $f1 | grep -i "0043 1083" | awk -F// '{print $3}' | cut -d'\' -f1`

if ($asset == "") then
echo  "no asset"
else
echo "asset factor = $asset"
set dwell = `echo "scale = 2; $dwell*$asset " | bc -l `
echo "Adjusted dwell time = $dwell usec"
endif

set dratio = `echo "scale = 10; $dwell/$tediff " | bc -l `
echo " "
echo "INFO: dratio = $dratio "


#Changed dicom2 to dicom_hdr
set l1 = `dicom_hdr $f1 | grep "0018 0080"`
set TR  = `echo "$l1" | awk -F// '{print $3}'`
echo ""
echo "INFO: EPI TR = $TR msec"
echo ""


if ($usrnifti) then
   echo "Use user-supplied NIFTI for EPI VOL"   
   cp  $epinifti $curdir/$tmpdir/epivol$postfix
   cd $curdir

else

echo "Creating NIFTI for EPI DICOM"
#note: use the find command since ls i* may not work when the number
#      of DICOM files is to large (depending on system defaults)

cd $curdir/$epidir

#KL added "sort" since "find" arranges file by the time stamp,
#and the time stamp is often reset duing data transfer. 
find . -type f -name 'i*' | sort > $curdir/$tmpdir/flist
set f1 = `cat $curdir/$tmpdir/flist | wc -l`
set nreps =  ` echo " $f1 / $nslices " | bc -l`

echo $f1 $nreps $nslices

to3d -prefix epivol$postfix -time:zt $nslices $nreps $TR alt+z -@ < $curdir/$tmpdir/flist

mv epivol$postfix $curdir/$tmpdir
cd $curdir

endif

# Check matrix size and resample fmap data if needed (KL):
set matfm = `fslinfo $tmpdir/$mag_te1 | awk '{if($1 == "dim1") print $2}; {if($1 == "dim2") print $2}; {if($1 == "dim3") print $2}'`


set matepi = `fslinfo $tmpdir/epivol$postfix | awk '{if($1 == "dim1") print $2}; {if($1 == "dim2") print $2}; {if($1 == "dim3") print $2}'`


if (($matfm[1] != $matepi[1]) || ($matfm[2] != $matepi[2]) || ($matfm[3] != $matepi[3])) then
echo "**************************** WARNING ****************************************"
echo "WARNING: fieldmap data ($matfm) does not match the functional data ($matepi) !"
echo "Resampling the fieldmaps to ($matepi) ... "
echo "*****************************************************************************"

set cmd=(3dresample -master  $tmpdir/epivol$postfix -prefix $tmpdir/mag_te1_rs$postfix -inset $tmpdir/$mag_te1)
echo $cmd
$cmd

set gp = $tmpdir/mag_te1_rs$postfix 

set cmd=(3dresample -master  $gp -prefix $tmpdir/re_dph_rs$postfix -inset $tmpdir/re_dph$postfix)
echo $cmd
$cmd

set cmd=(3dresample -master  $gp -prefix $tmpdir/im_dph_rs$postfix -inset $tmpdir/im_dph$postfix)
echo $cmd
$cmd

set cp_dph = $tmpdir/cp_dph_rs$postfix
fslcomplex -complex  $tmpdir/re_dph_rs$postfix $tmpdir/im_dph_rs$postfix $cp_dph

else
echo "check ... fieldmap data ($matfm) matches the functional data ... OK"
endif

echo "ppge done"
echo "calling epidewarp4.ucsd"


# call epidewarp now
    set cmd = (epidewarp4.ucsd --mag $gp --dph $cp_dph --epi $tmpdir/epivol$postfix \
               --tediff $tediff  --esp $dwell  --vsm $tmpdir/vsm$postfix --exfdw $tmpdir/ex$postfix \
               --epidw $outstem$postfix --nocleanup --postfix $postfix --nphase $nphase --betf $betf)

if ($domask) then
      set cmd = "$cmd --outmask";
endif

if ($domoco == 0) then
     set cmd = "$cmd --nomoco";
endif

if($#refnum > 0 ) then
    set cmd = "$cmd --refnum $refnum"
endif

if($#unwarpdir > 0 ) then
    set cmd = "$cmd --unwarpdir $unwarpdir"
endif

if($#fmap > 0 ) then
    set cmd = "$cmd --fmap $fmap"
endif

echo $cmd

$cmd
  

if ($#fmap > 0 ) then
  fslcpgeom $gp $fmap$postfix 
endif

#fslmerge unwarped volumes together

    if($domask) then
    fslmerge -t $outstem$postfix `find $tmpdir/muavwepi*$postfix`
    else
    fslmerge -t $outstem$postfix `find $tmpdir/uavwepi*$postfix`
    endif
    #copy orignal epi header to the dewarped epi since fugue does 
    #not always preserve the head correctly -KL
    fslcpgeom $tmpdir/epivol$postfix $outstem$postfix




if($docleanup) then
  if ($domoco == 1) then 
    mv  $tmpdir/*par $outstem.par
  endif
  
  echo "Deleting files in $tmpdir"
  rm -rf $tmpdir
endif


#add Rx info to the dataset header KL
set f1 = `ls $epidir/*.1`
set ff=`which dicomrx`
if (-e "$ff") then
   dicomrx -brik $outstem$postfix $f1
endif

exit 0;


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

  set flag = $argv[1]; shift;

  switch($flag)

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


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

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

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

    case "-TR":
       echo "ERROR: flag $flag is obsolete!"
       exit 1

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

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

    
    case "-nomask":
      set domask = 0;
      breaksw

    case "-nobrik":
      set dobrik = 0;
      breaksw

    case "-nomoco":
      set domoco = 0;
      breaksw

    case "-nocleanup":
      set docleanup = 0;
      breaksw

    case "-nogzip":
      set postfix = '.nii';
      breaksw

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

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


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

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

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

end

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

############--------------##################
check_params:

  if($#fm == 0) then
    echo "ERROR: must specify a FM directory"
    exit 1;
  endif

  if ( ! -e $fm) then
    echo "ERROR: $fm does not exist!"
    exit 1;
  endif


  if($#epidir == 0) then
    echo "ERROR: must specify an input EPI directory"
    exit 1;
  endif


  if ( ! -e $epidir) then
    echo "ERROR: $epidir does not exist!"
    exit 1;
  endif

  if ( $#epinifti != 0 ) then
     if (! -e $epinifti) then
        echo "ERROR: $epinifti does not exist!"
        exit 1;
     endif
  endif


  if($#outstem == 0) then
    echo "ERROR: must specify an output stem"
    exit 1;
  endif


  if($ncl == 0) then
    echo "ERROR: number of coil cannot be zero"
    exit 1;
  endif

  if($nphase == 0) then
    echo "ERROR: number of phase splits cannot be zero"
    exit 1;
  endif

 #if($betf <= 0 || $betf >= 1) then
 #   echo "ERROR: Value of betf must be between 0 and 1"
 #   exit 1;
 #endif

  if($#tediff == 0) then
    echo "INFO: tediff set to auto"
    set tediff = 'auto'
  endif


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

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


usage_exit:
    echo "Name"
    echo "     ggfm  - preprocesses GE DICOM FILES for unwarping and then"
    echo "              calls epidewarp4.ucsd"
    echo ""
    echo "System requirements"
    echo "     AFNI - AFNI_2006 and up, 32bit version"
    echo "     FSL  - FSL3.2 and up, FSLOUTPUTTYPE = NIFTI_GZ (or NIFTI if use -nogzip option)"
    echo ""
    echo "Synopsis"
    echo "     ggfm -d <fm data dir> -i<EPI dir> -o <outstem> [<options>]"
    echo ""
    echo "Required Arguments"
    echo "     -d <fm data dir>"
    echo "     -i <input EPI data dir>"
    echo "     -o <outstem>"
    echo "     Environment variable FSLOUTPUTTYPE need to be set to NIFTI_GZ or NIFTI "
    echo "     (use -nogzip option if NIFTI)"
    echo ""
    echo "Optional Arguments"
    echo "     -n <epi.nii.gz>:    users supplied epivol data (must be nii.gz format) " 
    echo "     -tmpdir <tmpdir>:   temporary file directory; default: <current directory>/tmp "
    echo "     -tediff <tediff>:   TE difference [us]; default: auto"
    echo "     -unwarpdir <dir>:   unwarping direction = x / y / z / x- / y- / z-, default = y"
    echo "     -avwvol         :   enables AVW output volume creation"
    echo "     -fmap           :   enables field map output volume creation"
    echo "     -nomask         :   disables brain masking of output EPI volume" 
    echo "     -nomoco         :   disables motion correction of EPI prior to unwarping"
    echo "     -refnum <refnum>:   Reference image number [0,nframes-1] for motion correction, default is middle image in a series"
    echo "     -nocleanup      :   disables removal of temporary files"
    echo "     -ncoil          :   number of receivers used to acquire fieldmap (needed only when using coils other than BODY and 8HRBRAIN)" 
    echo "     -nphase         :   number of phase split in phase unwrapping (default =8)" 
    echo "     -betf         :   fractional intensity threshold for BET (0->1) (default =0.5)" 
    echo "     -nogzip         :   don't gzip NIFTI files" 
    echo "" 
    echo "Outputs"
    echo "     <outstem> - unwarped volume filename stem" 
    echo ""
    echo "Version"
    echo "     "$VERSION
    echo ""
    echo "Credits"
    echo "     FSL library" 
    echo ""
    echo "Reporting Bugs"
    echo "     Report bugs to kunlu@ucsd.edu"
    echo ""



  if($PrintHelp) \
  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