#!/bin/csh -f # epidewarp2.ucsd # # This script uses FSL prelude and fugue tools to unwarp EPI images. # Original version is called epidewarp.fsl authored by Doug Greve at MGH for use # by the fBIRN consortium. # # This version is called epidewarp2.ucsd and is adapted from version 1.1 (2004/7/27) of # epidewarp.fsl. The intent is to faciliate unwarping from GE DICOM images at UCSD. # As much as possible, this script will be updated to reflect improvements # in epidewarp.fsl as these are made available. # # For processing at UCSD, this script is intended to be called by the script: ppge2 # # Version History # 1.0 TTL 040728 Added gemode input to support different processing for GE DICOM files at UCSD # TTL 040926 Checked into CVS repository # TTL 040927 Modifying to support multiple EPI frames # 3.0 TTL 040929 Adding motion correction features # 3.1 TTL 041001 Redefine refnum to go from 0 to nframes-1, in accordance with standard AVW indexing # 3.2 TTL 041014 Added unwarpdir option # 3.3 GTB 050413 Added write fieldmap option # 3.4 GTB 050505 Added -plots flag on mcflirt # 3.5 GTB 061505 Modified for NIFTI format # Send Comments/Questions to giedrius@salk.edu OR ttliu@ucsd.edu set VERSION = '$Id: epidewarp2.ucsd,v 1.2 2005/08/11 20:41:52 webtt Exp $' set inputargs = ($argv); set do_outmask = 0; set gemode = 0; set domoco = 1; set dph = (); set dph2 = (); set mag = (); set epi = (); set refnum = (); set unwarpdir = (); set fmap = (); # Difference between first and second echoes of the B0 Map set tediff = (); # Suggest 2.44 ms for 3T # EPI Echo Spacing (MGH .58 ms) set esp = (); # FUGUE parameters set sigma = (); set npart = (); # Outputs set epidw = (); set vsm = (); set exfdw = (); 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 exf = $tmpdir/exf.nii.gz set brain = $tmpdir/brain.nii.gz set head = $tmpdir/head.nii.gz if($epimerged) then # Extract the middle time point for the example func (exf) set nframes = `avwinfo $epi | awk '{if($1 == "dim4") print $2}'` if($nframes == 1) then set nmidframe = 0; else set nmidframe = `echo "$nframes/2-1" | 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; # See if there's a .mat file to propogate set epibase = `basename $epi .nii.gz`; set epimat = $epibase.mat if(! -e $epimat) set epimat = (); else set nframes = `ls -1 $epi*.nii.gz | wc -l` echo "INFO: $nframes EPI time points" set nmidframe = `echo "$nframes/2" | bc ` set exf = ` printf %s%04d.nii.gz $epi $nmidframe ` echo $exf # fix this later set epimat = () endif #REDEFINE exf if refnum exists if($#refnum > 0 ) then if ($refnum >= 0 && $refnum < $nframes) then if($epimerged) then set cmd = (avwroi $epi $exf $refnum 1) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; set nmidframe = $refnum else @ refnum = $refnum + 1 set exf = ` printf %s%04d.nii.gz $epi $refnum ` echo $exf endif else echo "refnum should be >= 0 and < $nframes" exit 1 endif 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.nii.gz # 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; # 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 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # MAKE THIS PART DIFFERENT FOR GE DICOM #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if ($gemode == 1) then #NOTE: In gemode the field map inputs are complex AVW volumes # with one volume per TE. # The phases for TE1 and TE2 data are put into ph1 and dph, respectively # These are later merged into ph2 # Do the phase unwrapping of first echo (-f for 3D, -v for verbose) set cmd = (prelude -c $dph -o $tmpdir/ph1.nii.gz -f -v -m $head); echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; set cmd = (prelude -c $dph2 -o $tmpdir/dph.nii.gz -f -v -m $head); echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; set dph = $tmpdir/dph.nii.gz set ph1 = $tmpdir/ph1.nii.gz else # 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 # 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.nii.gz set cmd = (avwmaths $dph -mul 0 $ph1) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; endif # Merge, baby, merge set ph2 = $tmpdir/ph2.nii.gz 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. set vsmmag = $tmpdir/vsmmag.nii.gz set magdw = $tmpdir/magdw.nii.gz # To be thrown away set cmd = (fugue -i $mag -u $magdw -p $ph2 \ --dwell=$esp --asym=$tediff --mask=$brain --saveshift=$vsmmag); if($#unwarpdir > 0 ) then set cmd = "$cmd --unwarpdir=$unwarpdir"; endif if($#fmap > 0 ) then set cmd = "$cmd --savefmap=$fmap"; endif 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.nii.gz set cmd = (fugue -i $mag -w $magfw --loadshift=$vsmmag --mask=$brain ) if($#unwarpdir > 0 ) then set cmd = "$cmd --unwarpdir=$unwarpdir"; endif echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; # Register magfw to example func. There are some parameters here # that may need to be tweeked. set cmd = (flirt -in $magfw -ref $exf \ -out $tmpdir/magfw-in-exf.nii.gz \ -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; # 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 .nii.gz`; cp $epimat $vsmbase.mat endif # Check whether we can stop at this point if($#exfdw == 0 && $#epidw == 0) goto done; # UNWARP ONLY THE EXAMPLE FRAME HERE if($#exfdw != 0) then # Now apply the VSM to the exf (um=unmasked) set exfdwum = $tmpdir/exfdwum.nii.gz set cmd = (fugue -i $exf -u $exfdwum --loadshift=$vsm --mask=$brain ); if($#unwarpdir > 0 ) then set cmd = "$cmd --unwarpdir=$unwarpdir"; endif echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; # Now mask the dewarped exf using the mag brain mask. # This is necessary to prevent voxels in the original EPI # that are out of the brain from simply being copied into # the dewarped image. set cmd = (avwmaths $exfdwum -mul $brain $exfdw) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; if($#epimat) then # Propogate mat file set exfbase = `basename $exfdw .nii.gz`; cp $epimat $exfbase.mat endif endif # UNWARP ALL EPI FRAMES HERE if($#epidw != 0) then if($domoco == 1 && $epimerged == 1) then set mepi = $tmpdir/mepi.nii.gz set cmd = (mcflirt -in $epi -out $mepi -refvol $nmidframe -plots ) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; set epi = $mepi endif @ frame = 1; while($frame <= $nframes) if ($epimerged) then set inimg = $tmpdir/inframe.nii.gz @ roiframe = $frame - 1 set cmd = (avwroi $epi $inimg $roiframe 1) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; else set inimg = `printf %s%04d.nii.gz $epi $frame`; if ($domoco) then set cmd = (flirt -in $inimg -ref $exf -out $tmpdir/moco.ni.gz -dof 6 \ -searchrx -10 10 \ -searchry -10 10 \ -searchrz -10 10 ) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; set inimg = $tmpdir/moco.nii.gz endif endif set outimg = `printf $tmpdir/uavwepi%04d.nii.gz $frame`; echo $inimg $outimg set cmd = (fugue -i $inimg -u $outimg --loadshift=$vsm --mask=$brain ); if($#unwarpdir > 0 ) then set cmd = "$cmd --unwarpdir=$unwarpdir"; endif echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; #OPTIONAL MASKING OF THE DATA if($do_outmask) then set moutimg = `printf $tmpdir/muavwepi%04d.nii.gz $frame`; set cmd = (avwmaths $outimg -mul $brain $moutimg) echo $cmd |& tee -a $LF $cmd |& tee -a $LF if($status) exit 1; endif # Need to add propagation of mat file! @ frame = $frame + 1; end 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.ucsd 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 "--dph2": if ( $#argv == 0) goto arg1err; set dph2 = $argv[1]; shift; if(! -e $dph2) then echo "ERROR: cannot find $dph2" exit 1; endif set gemode = 1 breaksw case "--epi": if ( $#argv == 0) goto arg1err; set epi = $argv[1]; shift; set epimerged = 1 if(! -e $epi) then echo "ERROR: cannot find $epi" exit 1; endif breaksw case "--sepi": if ( $#argv == 0) goto arg1err; set epi = $argv[1]; shift; set epimerged = 0 if(! -e {$epi}0001.nii.gz) then echo "ERROR: cannot find {$epi}0001.nii.gz (first image in series)" 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 "--fmap": if ( $#argv == 0) goto arg1err; set fmap = $argv[1]; shift; breaksw case "--unwarpdir": if ( $#argv == 0) goto arg1err; set unwarpdir = $argv[1]; shift; breaksw case "--tmpdir": if ( $#argv == 0) goto arg1err; set tmpdir = $argv[1]; shift; set cleanup = 0; breaksw case "--refnum": set refnum = $argv[1]; shift; breaksw case "--nomoco": set domoco = 0; breaksw case "--nocleanup": set cleanup = 0; breaksw case "--outmask": set do_outmask = 1; 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($#dph2 == 0) then echo "ERROR: no phase2 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($#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 mkdir -p $tmpdir if($status) then echo "ERROR: cannot create $tmpdir" exit 1; endif if($#epidw != 0) then 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: epidewarp2.ucsd" echo "" echo "Inputs" echo " --mag volid : B0 magnitude volume" echo " --dph volid : B0 phase difference volume OR complex data from echo 1 if gemode = 1" echo " --dph2 volid : complex data from echo 2; specifiying this input will set gemode = 1" echo " --epi volid : epi volume (note: specify this for merged AVW data)" echo " --sepi volid : epi series prefix (note: specify this for non-merged AVW data) " echo " --tediff tediff : difference in B0 field map TEs" echo " --esp esp : EPI echo spacing" echo " " echo "Outputs" echo " --vsm volid : voxel shift map (required)" echo " --exfdw volid : dewarped example func volume" echo " --epidw volid : dewarped epi volume" echo " --fmap volid : fieldmap volume" echo " " echo "Options " echo " --unwarpdir