:orphan: .. _ahelp_epi_b0_correct.py: ***************** epi_b0_correct.py ***************** .. contents:: :local: | PURPOSE ======= .. code-block:: none 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 ====== .. code-block:: none + 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 ======= .. code-block:: none + 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 ======= .. code-block:: none -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 applie d; 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 +++++++++++++++++++++++++++++++++ .. code-block:: none 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 ~~~~~~~~~~~~~~~~~~~~~ .. code-block:: none 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 ~~~~~~~~~~~~~~~~~~~~~~~ .. code-block:: none 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 ++++++++++++++++ .. code-block:: none 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 ======== .. code-block:: none # 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