#!/bin/csh -f
# epidewarp.fsl
#

set VERSION = '$Id: epidewarp.fsl,v 1.17 2005/06/08 00:55:43 greve Exp $'
set inputargs = ($argv);

set dph = ();
set mag = ();
set epi = ();

# Difference between first and second echoes of the B0 Map
set tediff = (); # Suggest 2.46 ms for 3T

# EPI Echo Spacing (MGH .58 ms)
set esp    = (); 

# FUGUE parameters
set sigmamm = 2; # Smoothing sigma in mm
set npart = ();

# Outputs
set epidw = ();
set vsm = ();
set vsmmag = ();
set exfdw = ();
set exf = ();
set DoMagExfReg = 1;

set tmpdir = ();
set cleanup = 1;
set cleanup_forced = 0;
set PrintHelp = 0;

## If there are no arguments, just print useage and exit ##
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

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

# Create a log file
set LF = $vsm.log
if(-e $LF) mv $LF $LF.bak
echo Logfile is $LF
date >> $LF
echo $VERSION >> $LF
pwd >> $LF
echo $0  >> $LF
echo $inputargs >> $LF
df -h $tmpdir  >> $LF
df -h $outdir  >> $LF
which prelude >> $LF
which fugue   >> $LF

# Some temp files
set brain = $tmpdir/brain.img
set head  = $tmpdir/head.img

# See if there's a .mat file to propogate
if($#epi) then
  set epibase = `basename $epi .img`;
  set epidir  = `dirname $epi`;
else
  set epibase = `basename $exf .img`;
  set epidir  = `dirname $exf`;
endif
set epimat = $epidir/$epibase.mat
if(! -e $epimat) set epimat = ();

if($#exf == 0) then
  # Extract the middle time point for the example func (exf)
  set exf   = $tmpdir/exf.img
  set nframes = `avwinfo $epi | awk '{if($1 == "dim4") print $2}'`
  if($nframes == 1) then
    set nmidframe = 0;
  else
    set nmidframe = `echo "$nframes/2" | bc `;
  endif
  echo "nframes = $nframes, nmidframe = $nmidframe" |& tee -a $LF
  set cmd = (avwroi $epi $exf $nmidframe 1)
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  if($#epimat) cp $epimat `dirname $exf`/`basename $exf .img`.mat 
endif
if(! -e $exf) then
  echo "ERROR: cannot find $exf"
  exit 1;
endif

# Keep only the first frame from mag
set cmd = (avwroi $mag $tmpdir/mag 0 1)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
set mag = $tmpdir/mag.img
if($#epimat) cp $epimat `dirname $mag`/`basename $mag .img`.mat 

# Create brain mask from the mag
set cmd = (bet $mag $brain)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
set cmd = (avwmaths $brain -bin $brain)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($#epimat) cp $epimat `dirname $brain`/`basename $brain .img`.mat 

# Create head mask by dilating the brain mask 3 times
@ nthdil = 1;
while($nthdil <= 3)
  if($nthdil == 1) then
    set cmd = (avwmaths $brain -dil $head)
  else
    set cmd = (avwmaths $head  -dil $head)
  endif
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  @ nthdil = $nthdil + 1;
end
if($#epimat) cp $epimat `dirname $head`/`basename $head .img`.mat # Copy matfile

# Rescale the delta phase to be between -pi and pi. Starts out
# at 0 - 4095. Make sure the phase is float precision with _32R
set cmd = (avwmaths_32R $dph -sub 2047.5 -mul 0.00153436 $tmpdir/dph)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
set dph = $tmpdir/dph
if($#epimat) cp $epimat `dirname $dph`/`basename $dph .img`.mat # Copy matfile

# Do the phase unwrapping (-f for 3D, -v for verbose)
set cmd = (prelude -p $dph -a $mag -o $dph -f -v -m $head);
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;

# FUGUE wants a phase image for each echo, but we only have
# phase difference between echoes. So create an image of 0s 
# and merging with the phase diff.
set ph1 = $tmpdir/ph1.img
set cmd = (avwmaths $dph -mul 0 $ph1)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
# Merge, baby, merge
set ph2 = $tmpdir/ph.img
set cmd = (avwmerge -t $ph2 $ph1 $dph)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;

# Create the voxel shift map (VSM) in the mag/phase space. Use mag as 
# input to assure that VSM is same dimension as mag. The input only affects
# the output dimension. The content of the input has no effect
# on the VSM. The dewarped mag volume is meaningless and will be thrown away. 
if(! $#vsmmag) set vsmmag = $tmpdir/vsmmag.img
set magdw = $tmpdir/magdw.img # To be thrown away
set cmd = (fugue -i $mag -u $magdw -p $ph2 \
    --dwell=$esp --asym=$tediff --mask=$brain --saveshift=$vsmmag);
if($#sigmamm && $sigmamm > 0) set cmd = ($cmd --smooth2=$sigmamm);
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($#epimat) cp $epimat `dirname $vsmmag`/`basename $vsmmag .img`.mat 

# Remove the mean in-brain shift from VSM
set vsmmean = `avwstats $vsmmag -M`;
set cmd = (avwmaths $vsmmag -sub $vsmmean $vsmmag)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
set cmd = (avwmaths $vsmmag -mul $brain $vsmmag);
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;

# Forward warp the mag in order to reg with func
# What does mask do here?
set magfw = $tmpdir/magfw.img
set cmd = (fugue -i $mag -w $magfw --loadshift=$vsmmag --mask=$brain )
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($#epimat) cp $epimat `dirname $magfw`/`basename $magfw .img`.mat

# Register magfw to example func. There are some parameters here
# that may need to be tweeked. Should probably strip the mag.
if($DoMagExfReg) then 
  set cmd = (flirt -in $magfw -ref $exf \
    -out $tmpdir/magfw-in-exf.img \
    -omat $tmpdir/magfw-in-exf.fsl.mat \
    -bins 256 -cost corratio \
    -searchrx -10 10 \
    -searchry -10 10 \
    -searchrz -10 10 \
    -dof 6  -interp trilinear)
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  if($#epimat) cp $epimat $tmpdir/magfw-in-exf.mat # Copy matfile
else
  echo "Not performing mag-exf registration "|& tee -a $LF
  echo "1 0 0 0" >  $tmpdir/magfw-in-exf.fsl.mat
  echo "0 1 0 0" >> $tmpdir/magfw-in-exf.fsl.mat
  echo "0 0 1 0" >> $tmpdir/magfw-in-exf.fsl.mat
  echo "0 0 0 1" >> $tmpdir/magfw-in-exf.fsl.mat
endif


# Now resample VSM into epi space. This will take care of any 
# differences in in-plane voxel size.
set cmd = (flirt -in $vsmmag -ref $exf -out $vsm \
  -init $tmpdir/magfw-in-exf.fsl.mat -applyxfm)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($#epimat) then # Propogate mat file
  set vsmbase = `basename $vsm .img`;
  cp $epimat `dirname $vsm`/$vsmbase.mat
endif

# Now resample brain mask into epi space. This will take care of any 
# differences in in-plane voxel size.
set brainexf = $tmpdir/brainexf.img
set cmd = (flirt -in $brain -ref $exf -out $brainexf \
  -init $tmpdir/magfw-in-exf.fsl.mat -applyxfm)
echo $cmd |& tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($#epimat) then # Propogate mat file
  set base = `basename $brainexf .img`;
  cp $epimat `dirname $brainexf`/$base.mat
endif

# Check whether we can stop at this point
if($#exfdw == 0 && $#epidw == 0) goto done;

#------------------ Dewarp the exf --------------------------#
if($#exfdw != 0) then
  # Now apply the VSM to the exf
  set cmd = (fugue -i $exf -u $exfdw --loadshift=$vsm --mask=$brainexf );
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;

  if($#epimat) then # Propogate mat file
    set exfbase = `basename $exfdw .img`;
    cp $epimat `dirname $exfdw`/$exfbase.mat
  endif
endif

#------------------ Dewarp the epis --------------------------#
if($#epidw != 0) then
  set nframes = `avwinfo $epi | awk '{if($1 == "dim4") print $2}'`;

  # Split the epi into multiple volumes
  set epipath = `pwd`/$epi
  mkdir -p $tmpdir/epi-split
  pushd $tmpdir/epi-split
  set cmd = (avwsplit $epipath)
  echo $cmd 
  $cmd 
  if($status) exit 1;
  popd

  # Go through each frame
  @ nthframe = 0
  while($nthframe < $nframes)
    set invol  = $tmpdir/epi-split/`printf vol%04d.img $nthframe`;
    set outvol = $tmpdir/epi-split/`printf voldw%04d.img $nthframe`;
    set cmd = (fugue -i $invol -u $outvol --loadshift=$vsm --mask=$brain );
    echo $cmd |& tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
    @ nthframe = $nthframe + 1;
  end

  # Merge them back - could they make this more convoluted???
  set cmd = (avwmerge -t $epidw `ls $tmpdir/epi-split/voldw*img`);
  echo $cmd |& tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;

  if($#epimat) then # Propogate mat file
    set epidwbase = `basename $epidw .img`;
    cp $epimat `dirname $epidw`/$epidwbase.mat
  endif

endif


goto done;

exit 0;
##############################################################

##############################################################
done:
  if($cleanup) then 
    echo "Deleting tmp dir $tmpdir" |& tee -a $LF
    rm -r $tmpdir
  endif
  date |& tee -a $LF
  echo "epidewarp.fsl done" |& tee -a $LF
exit 0;

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

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "--mag":
      if ( $#argv == 0) goto arg1err;
      set mag = $argv[1]; shift;
      if(! -e $mag) then
        echo "ERROR: cannot find $mag"
        exit 1;
      endif
      breaksw

    case "--dph":
      if ( $#argv == 0) goto arg1err;
      set dph = $argv[1]; shift;
      if(! -e $dph) then
        echo "ERROR: cannot find $dph"
        exit 1;
      endif
      breaksw

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

    case "--epi":
      if ( $#argv == 0) goto arg1err;
      set epi = $argv[1]; shift;
      if(! -e $epi) then
        echo "ERROR: cannot find $epi"
        exit 1;
      endif
      breaksw

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

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

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

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

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

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

    case "--nomagexfreg":
      set DoMagExfReg = 0;
      breaksw

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

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

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

    case "--cleanup":
      set cleanup_forced = 1;
      breaksw

    case "--debug":
      set verbose = 1;
      set echo = 1; # turns on terminal echoing
      breaksw

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

end

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

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

  if($#mag == 0) then
    echo "ERROR: no magnitude volume specified"    
    exit 1;
  endif

  if($#dph == 0) then
    echo "ERROR: no phase diff volume specified"    
    exit 1;
  endif

  if($#tediff == 0) then
    echo "ERROR: no TE diff specified"    
    exit 1;
  endif

  if($#esp == 0) then
    echo "ERROR: no Echo Spacing specified"    
    exit 1;
  endif

  if($#epi == 0 && $exf == 0) then
    echo "ERROR: must specify either epi or exf."
    exit 1;
  endif

  if($#vsm == 0) then
    echo "ERROR: no output VSM specified"    
    exit 1;
  endif

  if($cleanup_forced) set cleanup = 1;

  set outdir = `dirname $vsm`;
  mkdir -p $outdir
  if($status) then
    echo "ERROR: cannot create $outdir"
    exit 1;
  endif

  if($#tmpdir == 0) set tmpdir = $outdir/tmp-epidewarp.$$.fsl
  mkdir -p $tmpdir
  if($status) then
    echo "ERROR: cannot create $tmpdir"
    exit 1;
  endif

  if($#epidw != 0) then
    if($#epi == 0) then
      echo "ERROR: need --epi with --epidw"
      exit 1;
    endif
    set epidwdir = `dirname $epidw`;
    mkdir -p $epidwdir
    if($status) then
      echo "ERROR: cannot create $epidwdir"
      exit 1;
    endif
  endif

  if($#exfdw != 0) then
    set exfdwdir = `dirname $exfdw`;
    mkdir -p $exfdwdir
    if($status) then
      echo "ERROR: cannot create $exfdwdir"
      exit 1;
    endif
  endif

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

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

############--------------##################
usage_exit:
  echo ""
  echo "USAGE: epidewarp.fsl"
  echo ""
  echo "Inputs"
  echo "  --mag volid     : B0 magnitude volume"
  echo "  --dph volid     : B0 phase difference volume"
  echo "  --exf volid     : example func volume (or use epi)"
  echo "  --epi volid     : epi volume to unwarp"
  echo "  --tediff tediff : difference in B0 field map TEs"
  echo "  --esp esp       : EPI echo spacing"
  echo "  --sigma sigmamm : 2D spatial gaussing smoothing stddev (def 2mm)"
  echo "  "
  echo "Outputs"
  echo "  --vsm volid   : voxel shift map (required)"
  echo "  --vsmmag volid: voxel shift map in mag space"
  echo "  --exfdw volid : dewarped example func volume"
  echo "  --epidw volid : dewarped epi volume"
  echo "  "
  echo "  --nomagexfreg : assume mag and exf are in register"
  echo "  --tmpdir dir  : save intermediate results here"
  echo "  --nocleanup   : do not delete tmpdir"
  echo "  --cleanup     : force deletion of tmpdir"
  echo "  --debug"
  echo "  --help"
  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

Help Outline:
  SUMMARY
  ALGORITHM OVERVIEW
  ARGUMENTS
  SUGGESTED FIELD MAP PARAMETERS
  NOTES
  SEE ALSO
  BUGS
  BUG REPORTING
  AUTHORS
  REFERENCES AND ACKNOWLEDGMENTS

SUMMARY

Front end for FSLs PRELUDE and FUGUE programs to correct for B0
distortion in functional EPI scans. The programs use a B0 fieldmap.
This is assumed to be two conventional GREs collected at two different
echo times (TEs). The field map should be acquired with the same slice
prescription, slice thickness, slice skip, voxel resolution, and
field-of-view as the EPI. This script uses only FSL commands (ie,
FreeSurfer not needed).

For the stock Siemens field map, two field maps are required, one to
get the phase and another to get the magnitude (this will change with
new version of the Siemens scanner software). The volume that is
stored as the phase is actually the phase difference and is scaled
between 0 to 4095 for -pi to pi. The magnitude volumes for both TEs
are saved, but only one is needed.

All volumes are assumed to be in analyze 4D format. All volumes should
be referred to as volid.img. If the EPI or EXF volume has a .mat file,
this will be propagated to the outputs.

ALGORITHM OVERVIEW

1.  Create a brain mask from the mag volume (BET)
2.  Create a head mask by dilating the brain mask 3 times
3.  Rescale the phase image to -pi to pi
4.  Unwrap the phase (PRELUDE)
5.  Create and smooth the voxel shift map (VSM) (FUGUE)
6.  Remove in-brain mean from VSM
7.  Forward warp the mag volume (FUGUE)
8.  Register the forward warped mag with the example func (FLIRT)
9.  Resample the VSM into EPI space (FLIRT)
10. Dewarp the EPI and/or Example Func (FUGUE).

ARGUMENTS

--mag magvolid

Magnitude volume from the B0 fieldmap. If more than one frame is
present, then only the first is used.

--dph phasediffvolid

Phase difference volume (Echo2-Echo1). These are assumed to be scaled
between 0 to 4095 for -pi to pi. Eg, dph.img

--exf examplevolid

Single volume to use as the example func. If not supplied, then the
middle time point from the EPIs will be used. This volume must have
enough anatomy to be able to register with the mag. It should also
have the exact same geometry as any EPIs. This option allows the user
to specify a functional for the example and then a statistical map as
the EPI.

--epi epivolid

EPI volume to dewarp. This could be a time series or a single frame
statistical map. If it is a stat map, then supply an example func
(--exf) to assure that there is some anatomy for registration to
the mag volume.

--tediff tediff

Difference in the echo times of the B0 field map in ms. This number is
set when acquiring the field map. This is often set to 2.46 ms at 3T,
which is the amount of time it takes for the fat and water to
rephase. The rephase time scales with field strength.

--esp echospacing

Time (ms) between the start of the readout of two successive lines in
k-space during the EPI acquisition.  This can also be thought of as
the time between rows.  This parameter is not available in the Siemens
DICOM header. However, it can be obtained from the console.  Open the
protocol for the functional scan (NOT the field map). Go to the
Sequence tab. Go to Part 1. The Echo Spacing is displayed on the
bottom right. It will also be on the protocol print out (select the
functional protocol in the Exam Explorer, right click, and select
Print...).  If one saves the raw kspace data, then the echo spacing
can be found next to the m_lEchoSpacing tag in the meas.asc and is in
usec. Note that the echo spacing is not the same as the time to read
out a line of k-space and it cannot be computed from the bandwidth
because neither of these methods takes into account the dead time
between lines when no data are being read out.

--sigma sigmamm 

2D Smoothing VSM with gaussian with stddev of sigmamm (in mm). Default
is 2 mm. This tends to have minimal impact but can clean things up
around the edge of the brain.

--vsm vsmvolid

Voxel shift map. This is a volume the same size as the EPI. The value
at each voxel indicates the amount the voxel has shifted in the phase
encode direction due to B0 distortion. This argument is required.

--exfdw exfdwvolid

This is exf volume or the middle time point from the EPI time series
(the "example functional") with B0 distortion removed. Mainly good for
checking how good the dewarping is without having to dewarp the entire
time series.

--epidw epidwvolid

This is the EPI time series with B0 distortion removed.

--tmpdir tmpdir

Location to put the directory for storing temporary files. By default,
this will be in a directory called tmp-epidewarp.fsl under the
directory to hold the VSM volume. When --tmpdir is used or --nocleanup
is specified, then this directory will not be deleted. Otherwise or
with --cleanup it will automatically be deleted.

--nocleanup

Do not delete the tmp dir. --tmpdir automatically implies --nocleanup.

--cleanup

Forces deleting of the tmp dir regardless of whether --tmpdir or 
--nocleanup have been specified.

--debug

Prints copious amounts to the screen.

SUGGESTED FIELD MAP PARAMETERS

These field map parameters are being used by the Biomedical
Informatics Research Network (www.nbirn.net). For Siemens, we
recommend using the stock sequence called gre_field_mapping.  The
DeltaTE might differ slightly on individual systems from the values
listed below depending upon the exact field strength (eg, 2.89T
instead of exactly 3T).

TR = 500ms (1000ms with phased-array)
TE1: Minimum possible
TE2 = TE1 + DeltaTE

DeltaTE = 4.92ms for 1.5T
Flip Angle: 65 for 1.5T
DeltaTE = 2.46ms for 3T 
Flip Angle: 55 for 3T 
DeltaTE = 1.845ms for 4T
Flip Angle: 55 for 4T 

FatSat: Off
Bandwidth: Maximum
Slice Prescription, FOV, Thickness, Skip, Pixel Resolution, Matrix Size: 
same as functional

NOTES

The output will not necessarily be registered with the mag volume. To
compare the mag with the dewarped output, first register the mag 
to the dewarped output.

SEE ALSO

http://www.fmrib.ox.ac.uk/fsl/fugue
http://www.ami.hut.fi/instructions/epi_distortion.html

BUGS

None (yet).

BUG REPORTING

This software is offered as is with no promise or implication that it
will be supported. However, if you have a problem or question, you can
send it to analysis-bugs@nmr.mgh.harvard.edu and we will get to as we
can. If you want your question to have a chance of being answered,
include all of the following: your command-line, text printed out by
the program, log file, and a description of the problem. Of course, if
you have made it this far into the documentation, you probably alrealy
know these things.

AUTHORS

Doug Greve, Dave Tuch, Tom Liu, and Bryon Mueller with generous help
from the FSL crew (www.fmrib.ox.ac.uk/fsl) and the Biomedical
Informatics Research Network (www.nbirn.net).

REFERENCES AND ACKNOWLEDGMENTS:

For prelude:

M. Jenkinson. A Fast, Automated, N-Dimensional Phase Unwrapping
Algorithm Magnetic Resonance in Medicine, 49(1), 193-197, 2003.

For fugue:

M. Jenkinson. Improved Unwarping of EPI Volumes using Regularised B0
Maps International Conference on Human Brain Mapping - HBM2001.

Jezzard, Peter, and Balaban, Robert S. Correction for Geometric
Distortion in Echo Planar Images from B0 Field Variations. Mag Res
Med, 1995, 34:65-73.

For epidewarp.fsl we ask that the following be placed in the 
acknowledgement section of you papers:

  "Some of the scripts involved in this work were developed by 
   FBIRN (www.nbirn.net)."