#!/bin/tcsh -f
# fs-topup - sources
if(-e $FREESURFER_HOME/sources.csh) then
  source $FREESURFER_HOME/sources.csh
endif

set VERSION = '$Id$';
set scriptname = `basename $0`
setenv IGNORE_TAG_RAS_XFORM

# The topup config depends on the number of slices. Eg, if odd, then
# use b02b0_1.cnf; if even then b02b0_2.cnf = b02b0.cnf, and topup
# will die if you give it the wrong one. However, if you don't give topup
# a config file at all, it will not die, though I'm not sure if it will
# do the right thing
set config = ()
set configodd  = /usr/pubsw/packages/fsl/6.0.1/src/topup/flirtsch/b02b0_1.cnf
set configeven = /usr/pubsw/packages/fsl/6.0.1/src/topup/flirtsch/b02b0_2.cnf
set UseConfig = 1 # 

set outdir = ();
set b0dc = ()
set subject = ();
set ForceUpdate = 0
set threads = 1
set vols = ()
set DoQA = 0
set ndilations = 1
set DoMask = 1
set params = ()
set MultiFrame = 1 # 1=compute mean, 2=take first frame
# Total readout time in sec. If TRT=1, then b0 map will be in voxels,
# otherwise in Hz. The choice of TRT will have two implications: (1)
# if the number is wrong, then the Hz value will be wrong. (2) if
# combining with FSL eddy, then it needs to be consistent. Note that
# FSL eddy typically fails when TRT=1. The choice itself will not 
# have an effect on the B0-corrected maps.
set TRT = 1 

set tmpdir = ();
set cleanup = 1;
set LF = ();

set inputargs = ($argv);
set PrintHelp = 0;
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:
setenv OMP_NUM_THREADS $threads

set StartTime = `date`;
set tSecStart = `date '+%s'`;
set year  = `date +%Y`
set month = `date +%m`
set day   = `date +%d`
set hour   = `date +%H`
set min    = `date +%M`

mkdir -p $outdir/log
if($#subject) mkdir -p $outdir/reg
pushd $outdir > /dev/null
set outdir = `pwd`;
popd > /dev/null

if($#tmpdir == 0) then
  if(-dw /scratch)   set tmpdir = /scratch/tmpdir.fs-topup.$$
  if(! -dw /scratch) set tmpdir = $outdir/tmpdir.fs-topup.$$
endif
#mkdir -p $tmpdir

# Set up log file
if($#LF == 0) set LF = $outdir/log/fs-topup.Y$year.M$month.D$day.H$hour.M$min.log
if($LF != /dev/null) rm -f $LF
echo "Log file for fs-topup" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
ls -l $0  | tee -a $LF
echo "" | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
uname -a  | tee -a $LF
echo "pid $$" | tee -a $LF
if($?PBS_JOBID) then
  echo "pbsjob $PBS_JOBID"  >> $LF
endif
if($?SLURM_JOB_ID) then
  echo SLURM_JOB_ID $SLURM_JOB_ID >> $LF
endif

if($UseConfig && $#config == 0) then
  # Get number of slices
  set tmpfile = /tmp/fs-topup.$$
  mri_info --o $tmpfile --dim $vols[1] 
  if($status) exit 1
  set dim = `cat $tmpfile`
  rm $tmpfile
  set nslices = $dim[3];
  set IsOdd = `python -c "print($nslices%2)"`
  echo "Dim: $dim  IsOdd $IsOdd" | tee -a $LF
  if($IsOdd) then
    set config = $configodd
  else
    set config = $configeven
  endif
endif

# Check for multiple frames. If so, compute the mean.  Should also do
# motion correction? Might not be needed as these scans are usually
# pretty quick.
set nframes = ()
set vols2 = ()
@ n = 0
foreach v ($vols)
  @ n = $n + 1
  set nframesfile = $outdir/log/nframes$n.dat
  set ud = `UpdateNeeded $nframesfile $v`
  if($ud || $ForceUpdate) then
    set cmd = (mri_info --nframes --o $nframesfile $v)
    echo $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) goto error_exit
  endif
  set nf = `cat $nframesfile`
  set nframes = ($nframes $nf)
  if($nf > 1 && $MultiFrame) then
    set vmn = $outdir/vol$n.nii.gz
    set ud = `UpdateNeeded $vmn $v`
    if($ud || $ForceUpdate) then
      if($MultiFrame == 1) set cmd = (mri_concat  $v --mean --o $vmn)
      if($MultiFrame == 2) set cmd = (mri_convert $v --frame 0  $vmn)
      echo $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) goto error_exit
    endif
    set vols2 = ($vols2 $vmn)
  else
    if($nf > 1 && $#params == 0) then
      echo "ERROR: multiframe data needs either --p or allowing of multiframemean"
      exit 1
    endif
    set vols2 = ($vols2 $v)
  endif
end
set vols = ($vols2)

# Concatenate the two directions (will be more if multiple frames allowed)
set both = $outdir/both.nii.gz
set ud = `UpdateNeeded $both $vols`
if($ud || $ForceUpdate) then
  set cmd = (mri_concat $vols --o $both)
  echo "\n" $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit
endif

# This indicates the phase encode direction. Below is for AP then
# PA. See also BBR command.  When running vol2vol or vol2surf or bbr
# on the 2nd volume have to change the PE dir. Note that when used
# with FSL eddy, this needs to be consistent.
if($#params == 0) then
  set params = $outdir/params.dat
  if(! -e $params) then
    echo "0  1 0 $TRT" > $params
    echo "0 -1 0 $TRT" >> $params
  endif
else
  cp $params $outdir/log
endif

# Determine if TRT=1
set TRTis1 = `echo "sqrt((1-$TRT)*(1-$TRT)) < .0001" | bc -l`

set b0dcnomask = $outdir/b0dc.nomask.nii.gz # vox if TRT=1, Hz otherwise
set bothb0dc = $outdir/both.b0dc.nii.gz
set ud = `UpdateNeeded $b0dcnomask $both`
if($ud || $ForceUpdate) then
  # Run FSL topup
  set cmd = (fs_time topup --imain=$both --datain=$params --out=$outdir \
   --fout=$b0dcnomask --iout=$bothb0dc)
  if($#config) set cmd = ($cmd --config=$config)
  if($threads > 1) set cmd = ($cmd --nthr=$threads) # FSL>=6.0.6, not sure it works
  echo "fsl-topup-command"  | tee -a $LF # makes it easier to find in log file
  echo "\n" $cmd | tee -a $LF
  fs_time $cmd |& tee -a $LF
  if($status) then
    if($threads > 1) then
      echo "If this error has to do with the --nthr topup option, then make sure you are using FSL>=6.0.6" | tee -a $LF
    endif
    goto error_exit
  endif
  # annoyingly creates ${outdir}_fieldcoef.nii.gz ${outdir}_movpar.txt
  # Move them to the output folder. Add "topup_" as a prefix so that
  # other FSL commands can take $outdir/topup as the topup argument.
  mv ${outdir}_fieldcoef.nii.gz $outdir/topup_fieldcoef.nii.gz 
  mv ${outdir}_movpar.txt $outdir/topup_movpar.txt
endif

# Compute the mean of the two (or more) corrected volumes
set meanb0dc = $outdir/mean.b0dc.nii.gz
set ud = `UpdateNeeded $meanb0dc $bothb0dc`
if($ud || $ForceUpdate) then
  set cmd = (mri_concat $bothb0dc --mean --o $meanb0dc)
  echo "\n" $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit
endif

if($DoMask) then
  # Create a mask in the corrected space
  set maskb0dc = $outdir/mask.b0dc.nii.gz
  set ud = `UpdateNeeded $maskb0dc $meanb0dc`
  if($ud || $ForceUpdate) then
    set cmd = (fs_time mri_synthstrip -i $meanb0dc -m $maskb0dc -t $threads)
    echo "\n" $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) goto error_exit
    if($ndilations > 0) then
      set cmd = (mri_binarize --i $maskb0dc --min 0.5 --dilate $ndilations --o $maskb0dc)
      echo "\n" $cmd | tee -a $LF
      $cmd | tee -a $LF
      if($status) goto error_exit
    endif
  endif
  # Now mask the B0map volume
  set b0dc = $outdir/b0dc.nii.gz
  set ud = `UpdateNeeded $b0dc $maskb0dc $b0dcnomask`
  if($ud || $ForceUpdate) then
    set cmd = (mri_mask $b0dcnomask $maskb0dc $b0dc)
    echo "\n" $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) goto error_exit
  endif
else
  set b0dc = $b0dcnomask
endif

# Create vox and hz maps depending upon TRT
set vsm = $outdir/b0dc.vox.nii.gz
if($TRTis1) then
  # If TRT is 1.0, then b0dc is the voxel shift map (VSM)
  pushd $outdir
  ln -fs b0dc.nii.gz b0dc.vox.nii.gz
  popd
else
  # If TRT is not 1.0, then b0dc is in Hz
  pushd $outdir
  ln -fs b0dc.nii.gz b0dc.hz.nii.gz
  popd
  # Create a VSM by multiplying Hz by TRT.
  set ud = `UpdateNeeded $vsm $b0dc`
  if($ud || $ForceUpdate) then
    set cmd = (fscalc $b0dc -mul $TRT -o $vsm)
    echo "\n" $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) goto error_exit
  endif
endif

if($DoQA == 2) then
  # This is a comparison of how FS vs FSL apply the B0 map.  Compare
  # both.b0dc.nii (the fsl-corrected vols) to v1.qa and v2.qa.
  # both.b0dc.nii has two frames corresponding to v1 and v2.  Note:
  # vol2vol will not do jacobian correction whereas topup will, so you
  # might see some bright spots in the FSL output. When I examined
  # these more carefully, FS and FSL were both very close, but I
  # thought FS was better for v2 but FSL was better for v1; bbregister
  # mincost agrees. Not sure what is going on. The RMS diff in the
  # registrations is 0.4mm and 0.3mm, so not such a big deal
  # Not that useful so have to spec --qa2
  mkdir -p $outdir/qa
  foreach v (1 2) # might not be perfect with multiple frames
    set vb0dc = $outdir/qa/v$v.qa.nii.gz
    set ud = `UpdateNeeded $vb0dc $vols[$v] $b0dc` 
    if($ud || $ForceUpdate) then
      set cmd = (mri_vol2vol --mov $vols[$v] --targ $meanb0dc --regheader --vsm $vsm --o $vb0dc)
      if($v == 2) set cmd = ($cmd --vsm-pedir -2)
      echo "\n" $cmd | tee -a $LF
      $cmd |& tee -a $LF
      if($status) goto error_exit
    endif
  end
endif

if($#subject) then
  # If a subject was passed, then run BBR, mostly to do QA in that the mincost
  # should be lower after correction
  set mov = $meanb0dc
  set regmean = $outdir/reg/reg.mean.b0dc.lta
  set ud = `UpdateNeeded $regmean $mov` # should add b0dc
  if($ud || $ForceUpdate) then
    set cmd = (bbregister --s $subject --mov $mov --bold --reg $regmean --threads $threads)
    echo "\n" $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) goto error_exit
  endif
 
  if($DoQA) then
    foreach v (1 2) # might not be perfect with multiple frames
      foreach c (nob0dc b0dc)
        set reg = $outdir/reg/reg.v$v.$c.lta
        set mov = $vols[$v]
        set vsmopt = ""
        if($c == b0dc) set vsmopt = "--vsm $vsm"
        set ud = `UpdateNeeded $reg $mov` # should add b0dc
        if($ud || $ForceUpdate) then
          set cmd = (bbregister --s $subject --mov $mov --bold --reg $reg $vsmopt \
            --init-reg $regmean --threads $threads)
          if($v == 2) set cmd = ($cmd --vsm-pedir -2)
          echo "\n" $cmd | tee -a $LF
          $cmd |& tee -a $LF
          if($status) goto error_exit
        endif
      end # cor
    end # vol
  endif # QA
endif

#========================================================

# Cleanup
# if($cleanup) rm -rf $tmpdir

# Done
echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunMin = `echo $tSecRun/60|bc -l`
set tRunMin = `printf %5.2f $tRunMin`
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "Fs-Topup-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Fs-Topup-Run-Time-Min $tRunMin" |& tee -a $LF
echo "Fs-Topup-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "fs-topup Done" |& tee -a $LF
exit 0

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

############--------------##################
error_exit:
echo "ERROR:"

exit 1;
###############################################

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

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

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

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

    case "--sd":
      if($#argv < 1) goto arg1err;
      setenv SUBJECTS_DIR $argv[1]; shift;
      breaksw

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

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

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

    case "--i":
      if($#argv < 2) goto arg2err;
      set vols = ($argv[1] $argv[2]);
      shift;shift;
      breaksw

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

    case "--no-config":
      set UseConfig = 0
      breaksw

    case "--config-odd":
      set config = $configodd
      breaksw

    case "--config-even":
      set config = $configeven
      breaksw

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

    case "--multi-frame-mean":
      set MultiFrame = 1
      breaksw
    case "--multi-frame-first":
      set MultiFrame = 2
      breaksw
    case "--no-multi-frame":
      set DoMultiFrameMean = 0
      breaksw

    case "--qa":
      set DoQA = 1
      breaksw
    case "--qa2":
      set DoQA = 2
      breaksw

    case "--no-mask":
      set DoMask = 0
      breaksw

    case "--force":
     set ForceUpdate = 1
     breaksw
    case "--no-force":
     set ForceUpdate = 0
     breaksw

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

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      breaksw

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

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

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

    case "--debug":
      set verbose = 1;
      set echo = 1;
      breaksw

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

end

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

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

if($#outdir == 0) then
  echo "ERROR: must spec outdir"
  exit 1;
endif
if($#vols == 0) then
  echo "ERROR: must spec input volumes"
  exit 1;
endif
if($#subject) then
  if(! -e $SUBJECTS_DIR/$subject) then
    echo "ERROR: cannot find $subject"
    exit 1;
  endif
else
  if($DoQA) then
    echo "ERROR: must specify --s subject with --qa"
    exit 1;
  endif
endif

if($#config) then
  if(!-e $config) then
    echo "ERROR: cannot find $config"
    exit 1
  endif
  if($UseConfig == 0) then
    echo "ERROR: cannot spec both --config and --no-config"
    exit 1
  endif
endif

if($#params) then
  if(! -e $params) then
    echo "ERROR: cannot find $params"
    exit 1
  endif
endif

#if($DoQA && $TRT != 1) then
#  echo "ERROR: TRT = $TRT, cannot use --qa if TRT != 1"
#  exit 1
#endif

which topup >& /dev/null
if($status) then
  echo "ERROR: cannot find topup in the path."
  echo " Do you have FSL installed and in your path?"
  exit 1
endif

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

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

############--------------##################
usage_exit:
  echo ""
  echo "fs-topup - FS front end for FSL's topup (developed under FSL 6.0.5.1)"
  echo " --i vol1 vol2 : vol1 and vol2 have opposite phase encodings"
  echo " --o outdir"
  echo " --ndilations ndil : dilate mask"
  echo " --no-mask : do not mask the B0 map"
  echo " --s subject : register B0-corrected to subject (also needed for QA using mincost)"
  echo " --config topup.cfn : topup config file (does not use a config file by default)"
  echo " --config-{odd,even}: use b02b0_{1,2}.cnf"
  echo " --p params : parameter file passed to topup (default is assuming AP/PA)"
  echo " --no-multi-frame-mean : do not comute the mean of multi-frame inputs (will need a param file)"
  echo " --qa : bbregister with and without B0DC to assure that the cost goes down (needs --s)"
  echo  "  Note: this can add 5-10min onto this script's execution time"
  echo " --trt total readout time in sec (default is $TRT)" 
  echo " --multi-frame-mean : compute the mean of multi-frame inputs"
  echo " --multi-frame-first : extract the first frame of multi-frame inputs"
  echo " --no-multi-frame : do not do anything to multi-frame inputs (requires --p)"
  echo " --threads nthreads : does not appear to help much with topup"
  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

This program is a FreeSurfer frontend for for FSLs topup B0
distortion correction command (developed under FSL 6.0.5.1). It
collects the inputs, runs topup, creates a mask using synthstrip,
registers to FreeSurfer subject (if supplied). See also fs-eddytopper.

There will be (possibly) three output B0 maps: b0dc.nii.gz,
b0dc.vox.nii.gz, b0dc.hz.nii.gz.  If TRT is set to a value different
than 1, then b0dc.nii.gz will be in Hz and b0dc.hz.nii.gz will be a
symlink pointing to b0dc.nii.gz. Note that TRT must be set accurately
for the Hz values to be accurate (see below).  If TRT=1, then
b0dc.nii.gz will be in voxels and b0dc.vox.nii.gz will be a symlink
pointing to b0dc.nii.gz (and b0dc.hz.nii.gz) will not exit. Note that
b0dc.vox.nii.gz will always be right regardless of TRT. Note that
these maps will have been masked unless --no-mask is used. There is an
output b0dc.nomask.nii.gz (in either Hz or vox depending upon TRT).

Setting the total readout time (TRT). The proper value will be TRT =
EPIFactor * EchoSpacing * PartialFourier / AccelFactor, where the
EchoSpacing is in seconds.  Eg, TRT = 116*0.58ms*.001/2 = 0.0336
seconds. If the TRT is not set correctly, the only implication is that
the Hz values will be wrong.  It will still apply the proper
correction and the voxel shift map will still be right. Note that SMS
acceleration does not reduce the effective echo spacing.

Note that topup will register the second input volume (ie, the
reversed phase encode direction) to the first and the output maps will
be in the (undistored) space of the first volume. This allows for the
reversed phase encoded image to be acquired at a different time as the
series to be corrected. For this application, one would then take the
first frame from the series as the non-reversed input (first input
volume). 

If you are going to use the output in conjunction with other
FreeSurfer programs, eg, FSFAST or bbregister, then pass those
programs the b0dc.vox.nii.gz as the --vsm argument. Within the FSFAST
structure, copy or link b0dc.vox.nii.gz to sess/fsd/b0dcmap.nii.gz.

If you are going to use the output with other FSL programs, then 
pass outdir/topup as the topup input.

When using in conjunction with FSLs eddy for DTI analysis, the first
volume input to topup must be the same as the first volume in --imain
to eddy_correct and the TRT must be consistent (as does the
acquisition parameter file). According to the eddy FAQ
https://fsl.fmrib.ox.ac.uk/fsl/docs/diffusion/eddy/FAQ/index.html the
TRT value does not matter that much. As long as it is consistent
between topup and eddy, they give indistinguishable results. However,
eddy may crash if TRT is too low (~<.01) or too high (~>.9).  Also see
fs-eddytopper.

The topup config depends on the number of slices. Eg, if odd, then use
b02b0_1.cnf; if even then b02b0_2.cnf = b02b0.cnf, and topup will die
if you give it the wrong one. By default, fs-topup will look at the
number of slices and choose the right one. Note: if you dont give
topup a config file at all, it will not die, though Im not sure if it
will do the right thing.

When passing a subject, it will create outdir/reg/reg.mean.b0dc.lta
which is the BBR registration of the mean of the two B0-corrected
inputs.

Quality Assurance (QA): if you add the --qa option with --s, then
fs-topup will register each of the input volumes to the FreeSurfer
anatomical using bbregister with and without B0 correction. If the B0
correction is working, then the best BBR cost should decrease. There
will be a reg folder in the output folder with the followinng files:

reg.v1.b0dc.dat.mincost
reg.v1.nob0dc.dat.mincost
reg.v2.b0dc.dat.mincost
reg.v2.nob0dc.dat.mincost

v1 is for the first input volume and v2 for the second.  The first
value in each file is the minimum cost found by BBR. The value in the
b0dc file should be less than that in the corresponding nob0dc
file. If not, then something went wrong. This is something you must
check yourself; fs-topup will not throw an error. These registrations
may add 5-10min on the execution time for fs-topup, so you may want to
try it on a few cases to convince yourself it is doing the right thing
then leave it off.

If either vol1 or vol2 has multiple frames, then the mean across
frames is computed and used. The first frame can be used instead by
specifying --multi-frame-first. 

CITATION: From: https://fsl.fmrib.ox.ac.uk/fsl/docs/diffusion/topup/index.html

If you use topup in your research, please make sure that you reference
at least the first of the articles listed below, and ideally both.
For your convenience, we provide example text (short and more details
versions), which you are welcome to use in your methods description.

Brief summary text: "Data was collected with reversed phase-encode
blips, resulting in pairs of images with distortions going in opposite
directions. From these pairs the susceptibility-induced off-resonance
field was estimated using a method similar to that described in
[Andersson 2003] as implemented in FSL [Smith 2004] and the two images
were combined into a single corrected one."

[Andersson 2003] J.L.R. Andersson, S. Skare, J. Ashburner How to
correct susceptibility distortions in spin-echo echo-planar images:
application to diffusion tensor imaging. NeuroImage, 20(2):870-888,
2003.

[Smith 2004] S.M. Smith, M. Jenkinson, M.W. Woolrich, C.F. Beckmann,
T.E.J. Behrens, H. Johansen-Berg, P.R. Bannister, M. De Luca,
I. Drobnjak, D.E. Flitney, R. Niazy, J. Saunders, J. Vickers,
Y. Zhang, N. De Stefano, J.M. Brady, and P.M. Matthews. Advances in
functional and structural MR image analysis and implementation as
FSL. NeuroImage, 23(S1):208-219, 2004.