epi_b0_correct.py¶
PURPOSE¶
This program performs B0 distortion correction along the phase encode
(PE) direction, using an acquired frequency (phase) image. It was
initially written by Vinai Roopchansingh (NIMH, NIH).
Ver : 2.64
Date : Sep 23, 2021
INPUTS¶
+ frequency dset : (req) phase volume, which should be of similar
spatial resolution/FOV of EPI dset to which it
will be applied. Expected units are:
angular freq = rad/s = 2*PI*(freq in Hz).
If your dataset is in different units,
you can apply an appropriate scaling via the
command line, as discussed in the 'NOTES', below.
+ EPI dset : (req) EPI dset to which the B0 distortion correction
is applied.
+ mask dset : (req) binary mask of subject's brain
OR
+ magnitude dset : (req) volume in same space as frequency dset for
automasking, to create brain mask; also useful for
QC imaging (even if a mask is input separately)
+ PE parameters : (req) a number of parameters related to the
EPI vol are required to be input, such as its
- PE direction (AP, PA, RL, etc.)
- bandwidth per pixel OR effective TE
Optional scaling can be applied to the freq dset
(e.g., if units need to be adjusted appropriately).
These parameters can be provided either
individually, or by providing an accompanying JSON
that might/should contain all necessary
information.
NB: If you input a parameter on the command line,
it will take precedence over one found in the
EPI's JSON, if you are also using that. Thus, if
you know the JSON has *wrong* information, you
can selectively ignore that when running this
program.
OUTPUTS¶
+ WARP dset : a file called PREFIX_WARP.nii.gz, containing
the warp along the phase encode axis (on the
EPI dset's grid, with its obliquity info)
+ script of commands : a script of the commands used to generate the
WARP dset (and EPI)
+ text file of params : a text file of parameters either input or
derived from inputs and the dsets. This is
useful for verifying the consistency of
analysis (i.e., as a sanity check). Can be
converted to a JSON, if needed. Units are
given for all; the 'Warp (mm) in mask,
20-100 %ile' field might be the most cryptic
entrant-- it is a histogram of values of the
final warp field within the mask, at the
20th, 40th, 60th, 80th and 100th %iles.
Cryptic no more!
+ EPI (un)warped dset : the EPI dset with the estimated distortion
correction applied to it (and obliquity info
matching the original EPI's); hopefully it
is unwarped...
+ QC image dir : a directory called PREFIX_QC/, containing
some (hopefully) useful QC images of both the
corrected and uncorrected EPI on the magn dset,
as well as the mask on the magn dset. All images
are shown in the coordinates of the EPI, whether
the EPI is in oblique or scanner coordinates (the
other dsets will have been transformed or "sent"
to those coords).
RUNNING¶
-prefix PP : (req) prefix of output files; can include path
-in_freq DSET_FREQ : (req) phase dset (frequency volume). Should
be of similar spatial resolution and FOV as
EPI dset to which it will be applied; also,
must be scaled appropriately, where the
expected units are: Hz.
-in_epi DSET_EPI : (req) EPI dset to which the B0 distortion
correction that I have spent so much time
calculating will be applied
-in_mask DSET_MASK : (req) mask of brain volume
or
-in_magn DSET_MAGN : (req) magnitude dset from which to estimate brain
mask; it can be useful to enter a magn dset even
if a mask is input, in order to use it as a
reference underlay in the QC image directory
-in_anat DSET_ANAT : (opt) if input, this dset will be used to make
the underlay for the automatically generated
QC images; if this dset is not provided, then
the DSET_MAGN will be used (and if that is
not provided, then the QC images will just
have the given EPI(s) as ulay-only)
-in_epi_json FJSON : (opt) Several parameters about the EPI
dset must be known for processing; these MIGHT
be encoded in a JSON file accompanying the
EPI dset. If so, you can input the file
and let The Program try to find+interpret them.
At present, desirable keys/tags in the JSON
(with the keyword args you would otherwise use
when running this program) are:
PhaseEncodingDirection (or use '-epi_pe_dir')
and then either of the following:
BandwidthPerPixelPhaseEncode (or use '-epi_pe_bwpp')
OR
EffectiveEchoSpacing (or use '-epi_pe_echo_sp')
-epi_pe_dir DD : (req) direction (axis) of phase encoding,
e.g., AP, PA, RL, ...
NB: the order matters, providing the PE direction
(and not just PE axis); thus, 'AP' implies the
PE direction is A>>P, and 'PA' that it is P>>A, etc.
(Can come from EPI's JSON; see '-in_epi_json'.)
-epi_pe_bwpp BW : (req) bandwidth per pixel (in Hz) in the EPI
dset along the phase encode direction.
(Can come from EPI's JSON; see '-in_epi_json'.)
OR
-epi_pe_echo_sp ES : (req) *effective* TE spacing of phase encoded
volume, in units of 's'
(Can come from EPI's JSON; see '-in_epi_json'.)
-epi_pe_voxdim FOV : (opt) voxel size along the EPI dset's phase
encode axis, in units of 'mm'; should just be
determined internally from the EPI dataset
-scale_freq SF : (opt) scale to apply to frequency volume,
for example to change units to match.
NB: a negative value would invert the warp
(probably would not want that...?) See the
'NOTES ..' below for more information about
scaling, esp. for particular vendors or known
units, like physical frequency (Hz). (def: SF=1.0)
-out_cmds OC : (opt) name of output script, recording
commands that were run during the processing
(def: script is output to file using entered
prefix PP: PP_cmds.tcsh). If user uses
this option, then 'OC' is treated as the full
filename, including path
-out_pars OP : (opt) name of output parameters, recording
some relevant values that were input, found or
calculated during the processing; the file is
a colon-separated list that can be turned
into a JSON with abids_json_tool.py, if desired.
(def: pars are output to file using entered
prefix PP: PP_pars.txt). If user uses
this option, then 'OP' is treated as the full
filename, including path
-wdir_name WD : working directory name (no path, will be located
in directory with output dsets); if not
provided, will be given automatic name,
starting '__work_B0_corr_' and ending with a
random alphanumeric string, e.g.,
'__work_B0_corr__9huoXQ7c0AV'
-blur_sigma BS : amount of blurring to apply to masked, phase
encode dset (def: BS = 9)
-do_recenter_freq MC : method for 3dROIstats to recenter the phase
(=freq) volume within the brain mask.
If the value of MC is 'NONE', then the phase
dset will not be recentered.
If the value of MC is some number (e.g.,
60.704), then the phase dset will be
recentered by this specific value (must be in
units of the original, input phase dset).
If you want to recenter by the mean value,
then the value of MC should be "MEAN" (all
capital letters): this is because 3dROIstats
doesn't take a "-mean" option (it is actually
the default there), so one is entering a flag
to be interpreted, not a literal opt name.
(def: MC = mode; NB: this
method can't be used if the input dset type
is float, at which point the program will
exit and whine at the user to choose another
method, such as 'MEAN')
-mask_dilate MD1 MD2 ...
: if automasking a magnitude image to create a
brain mask, one can provide 3dmask_tool-style
erosion and dilation parameters for the mask.
NB: this ONLY applies if masking a magn image,
not if you have just put in a mask (you can
dilate that separately with 3dmask_tool).
Typically, one might input two values here, with
MD1 being negative (to erode) and MD2 being
positive (to dilate).
(def: MD* = -2 1)
-no_clean : don't remove the temporary directory of intermed
files
-qc_box_focus_ulay : an option about the QC image output-- this will
have @chauffeur_afni use the option+value:
'-box_focus_slices AMASK_FOCUS_ULAY'
which focuses the montage slices views on an
automask of the ulay dset involved (typically the
magn or anat dset; might not be desirable if
neither is used, because then the ulay will be
either uncorrected and corrected EPIs, which will
have slightly different automasks and therefore
slightly different slices might be shown, making
comparisons more difficult)
-no_qc_image : don't make pretty QC images (why not??)
-help : display program help in terminal (consider
'-hview' to open help in a separate text editor)
-ver : display program version number in terminal
-date : display date of program's last update in terminal
NOTES¶
Units of frequency/phase/fieldmap¶
It is important to have your input phase/frequency volume contain
the correct units for this program. Here, we expect them to be in
units of angular frequency: "radians/second" (rad/s).
Re. fieldmaps in Hz¶
If your frequency map has units of physical frequency, 'cycles per
second' (= Hz), then you just provide a command line argument to
internally scale your data to the appropriate angular frequency
unit we desire to use.
Physicists tell us that angular frequency 'w' is related to
physical frequency 'f' as follows:
w = 2 * PI * f
~ 6.2831853 * f
Therefore, if you are *sure* that your frequency (phase) volume is
really in units of Hz, then you can use the following command line
argument to set things right for using it here:
'-scale_freq 6.2831853'
Not too painful!
Re. Siemens fieldmaps¶
If your frequency map is one output by Siemens, then consider the
following (but doublecheck that it really applies to your darling
dataset!):
The standard range of fieldmap values in that case appears to be
either [-4096, 4095] or [0, 4095], depending on how your data were
converted. You can check the range on your dset with, e.g.:
3dinfo -dmin -dmax FREQ_DSET
will will likely *approximately* match one of those ranges.
These ranges come from dividing the measured phases by 2*PI (one
full phase) and then multiplying by either 2*4096 or 4096,
respectively. One could multiply by that inverse ratio, putting
the dataset into units of radians ('rad'); however, we ultimately
want the input frequency volume to be in units of angular
frequency: 'rad/s' ('2*PI*Hz'). Therefore, we also want to divide
by the frequency dset's echo time difference; this could be
calculated from 'EchoTime1' and 'EchoTime2' in the freq dset's
JSON sidecar (or possibly provided directly as
'EchoTimeDifference' there). For example, the standard value of
this at 3T is about 2.46 ms (= 0.00246 s), but check what it is in
your own data!
*Therefore*, in many cases of Siemens 3T data, one should be able
to convert the scaled freq dset into the the desired units of ang
freq by scaling the fieldmap by 2*PI/(2*4096*0.00246) ~ 0.311785
or by 2*PI/(4096*0.00246) ~ 0.623569, respectively. This could be
done using, say, 3dcalc to make a new freq dset; or, you could
provide this magic value to the present command with the scaling
option:
FREQ DSET ~RANGE (potential) PROGRAM OPTION
---------------- --------------------------
[-4096, 4095] : '-scale_freq 0.311785'
[0, 4095] : '-scale_freq 0.623569'
It is worth repeating: be sure that these numbers *really* apply to
your data!
Output QC images¶
QC images are automatically generated and put into a subdirectory
called PREFIX_QC/. Images are provided as montages in each of the
axi, sag and cor planes; data are shown in the EPI coords (oblique
if the EPI were oblique). The QC sets have the following simple
names (NB: if one inputs an anat vol via '-anat ..', then the
'anat' replaces 'magn' in the following lists-- even in the QC
image filenames):
Names if there is a magn vol included
qc_00_ee_magn+mask = ulay: edge-enhanced magn
olay: mask dset
qc_01_ee_magn+iepi = ulay: edge-enhanced magn
olay: input EPI[0] (uncorr)
qc_02_ee_magn+oepi = ulay: edge-enhanced magn
olay: output EPI[0] (corr)
Names if there is NOT a magn vol included
qc_11_iepi = ulay: input EPI[0] (uncorr)
qc_12_oepi = ulay: output EPI[0] (corr)
EXAMPLES¶
# Ex 1: With mask supplied, created earlier from magnitude image
epi_b0_correct.py \
-epi_pe_echo_sp 0.00031 \
-epi_pe_dir AP \
-in_freq sub-001_frequency.nii.gz \
-in_mask sub-001_magnitude_MASK.nii.gz \
-in_epi epiRest-sub-001.nii.gz \
-prefix b0_corr
# Ex 2: Input *magnitude* dset, from which to calculate mask
epi_b0_correct.py \
-epi_pe_echo_sp 0.00031 \
-epi_pe_dir AP \
-in_freq sub-001_frequency.nii.gz \
-in_magn sub-001_magnitude.nii.gz \
-in_epi epiRest-sub-001.nii.gz \
-prefix b0_corr
# Ex 3: Same as above, but freq dset was in units of Hz (convert
# to angular freq, scaling by 2*PI~6.283185)
epi_b0_correct.py \
-epi_pe_echo_sp 0.00031 \
-epi_pe_dir AP \
-scale_freq 6.283185 \
-in_freq sub-001_frequency.nii.gz \
-in_magn sub-001_magnitude.nii.gz \
-in_epi epiRest-sub-001.nii.gz \
-prefix b0_corr
# Ex 4: Input a JSON file (sidecar) accompanying the freq volume,
# and hope that it has all the necessary parameters/fields for
# this program.
epi_b0_correct.py \
-in_epi_json sub-001_frequency.json \
-in_freq sub-001_frequency.nii.gz \
-in_magn sub-001_magnitude.nii.gz \
-in_epi epiRest-sub-001.nii.gz \
-prefix b0_corr
# Ex 5: Same as Ex 4, but include the anatomical as an underlay
# in the QC imaging, and have the snapshot program focus just
# on an automask region of that anat volume
epi_b0_correct.py \
-in_epi_json sub-001_frequency.json \
-in_freq sub-001_frequency.nii.gz \
-in_magn sub-001_magnitude.nii.gz \
-in_epi epiRest-sub-001.nii.gz \
-in_anat sub-001_run-02_T1w+orig.HEAD \
-qc_box_focus_ulay \
-prefix b0_corr
