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