Susceptibility Distortion Correction (SDC)

Introduction

SDC methods usually try to make a good estimate of the field inhomogeneity map. The inhomogeneity map is directly related to the displacement of a given pixel \((x, y, z)\) along the PE direction (\(d_\text{PE}(x, y, z)\)) is proportional to the slice readout time (\(T_\text{ro}\)) and the field inhomogeneity (\(\Delta B_0(x, y, z)\)) as follows ([Jezzard1995], [Hutton2002]):

\[d_\text{PE}(x, y, z) = \gamma \Delta B_0(x, y, z) T_\text{ro} \qquad (1)\]

where \(\gamma\) is the gyromagnetic ratio. Therefore, the displacements map \(d_\text{PE}(x, y, z)\) can be estimated either via estimating the inhomogeneity map \(\Delta B_0(x, y, z)\) (Phase-difference B0 estimation and Direct B0 mapping sequences) or via image registration (Phase Encoding POLARity (PEPOLAR) techniques, Fieldmap-less estimation (experimental)).

Correction methods

The are five broad families of methodologies for mapping the field:

  1. Phase Encoding POLARity (PEPOLAR) techniques (also called blip-up/blip-down): acquire at least two images with varying PE directions. Hence, the realization of distortion is different between the different acquisitions. The displacements map \(d_\text{PE}(x, y, z)\) is estimated with an image registration process between the different PE acquisitions, regularized by the readout time \(T_\text{ro}\). Corresponds to 8.9.4 of BIDS.
  2. Direct B0 mapping sequences: some sequences (such as SE) are able to measure the fieldmap \(\Delta B_0(x, y, z)\) directly. Corresponds to section 8.9.3 of BIDS.
  3. Phase-difference B0 estimation: to estimate the fieldmap \(\Delta B_0(x, y, z)\), these methods measure the phase evolution in time between two close GRE acquisitions. Corresponds to the sections 8.9.1 and 8.9.2 of the BIDS specification.
  4. Fieldmap-less estimation (experimental): FMRIPREP now experimentally supports displacement field estimation in the absence of fieldmaps via nonlinear registration.
  5. Point-spread function acquisition: Not supported by FMRIPREP.

In order to select the appropriate estimation workflow, the input BIDS dataset is first queried to find the available field-mapping techniques (see Automatic selection of the appropriate SDC method). Once the field-map (or the corresponding displacement field) is estimated, the distortion can be accounted for (see Unwarping).

Calculating the effective echo-spacing and total-readout time

To solve (1), all methods (with the exception of the fieldmap-less approach) will require information about the in-plane speed of the EPI scheme used in acquisition by reading either the \(T_\text{ro}\) (total-readout time) or \(t_\text{ees}\) (effective echo-spacing):

fmriprep.interfaces.fmap.get_ees(in_meta, in_file=None)[source]

Calculate the effective echo spacing \(t_\text{ees}\) for an input EPI scan.

There are several procedures to calculate the effective echo spacing. The basic one is that an EffectiveEchoSpacing field is set in the JSON sidecar. The following examples use an 'epi.nii.gz' file-stub which has 90 pixels in the j-axis encoding direction.

>>> meta = {'EffectiveEchoSpacing': 0.00059,
...         'PhaseEncodingDirection': 'j-'}
>>> get_ees(meta)
0.00059

If the total readout time \(T_\text{ro}\) (TotalReadoutTime BIDS field) is provided, then the effective echo spacing can be calculated reading the number of voxels \(N_\text{PE}\) along the readout direction and the parallel acceleration factor of the EPI

\[= T_\text{ro} \, (N_\text{PE} / f_\text{acc} - 1)^{-1}\]

where \(N_y\) is the number of pixels along the phase-encoding direction \(y\), and \(f_\text{acc}\) is the parallel imaging acceleration factor (GRAPPA, ARC, etc.).

>>> meta = {'TotalReadoutTime': 0.02596,
...         'PhaseEncodingDirection': 'j-',
...         'ParallelReductionFactorInPlane': 2}
>>> get_ees(meta, in_file='epi.nii.gz')
0.00059

Some vendors, like Philips, store different parameter names (see http://dbic.dartmouth.edu/pipermail/mrusers/attachments/20141112/eb1d20e6/attachment.pdf):

>>> meta = {'WaterFatShift': 8.129,
...         'MagneticFieldStrength': 3,
...         'PhaseEncodingDirection': 'j-',
...         'ParallelReductionFactorInPlane': 2}
>>> get_ees(meta, in_file='epi.nii.gz')
0.00041602630141921826
fmriprep.interfaces.fmap.get_trt(in_meta, in_file=None)[source]

Calculate the total readout time for an input EPI scan.

There are several procedures to calculate the total readout time. The basic one is that a TotalReadoutTime field is set in the JSON sidecar. The following examples use an 'epi.nii.gz' file-stub which has 90 pixels in the j-axis encoding direction.

>>> meta = {'TotalReadoutTime': 0.02596}
>>> get_trt(meta)
0.02596

If the effective echo spacing \(t_\text{ees}\) (EffectiveEchoSpacing BIDS field) is provided, then the total readout time can be calculated reading the number of voxels along the readout direction \(T_\text{ro}\) and the parallel acceleration factor of the EPI \(f_\text{acc}\).

\[T_\text{ro} = t_\text{ees} \, (N_\text{PE} / f_\text{acc} - 1)\]
>>> meta = {'EffectiveEchoSpacing': 0.00059,
...         'PhaseEncodingDirection': 'j-',
...         'ParallelReductionFactorInPlane': 2}
>>> get_trt(meta, in_file='epi.nii.gz')
0.02596

Some vendors, like Philips, store different parameter names:

>>> meta = {'WaterFatShift': 8.129,
...         'MagneticFieldStrength': 3,
...         'PhaseEncodingDirection': 'j-',
...         'ParallelReductionFactorInPlane': 2}
>>> get_trt(meta, in_file='epi.nii.gz')
0.018721183563864822

From the phase-difference map to a field map

To solve (1) using a phase-difference map, the field map \(\Delta B_0(x, y, z)\) can be derived from the phase-difference map:

fmriprep.interfaces.fmap.phdiff2fmap(in_file, delta_te, newpath=None)[source]

Converts the input phase-difference map into a fieldmap in Hz, using the eq. (1) of [Hutton2002]:

\[\Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}}\]

In this case, we do not take into account the gyromagnetic ratio of the proton (\(\gamma\)), since it will be applied inside TOPUP:

\[\Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}}\]

References

[Jezzard1995]P. Jezzard, R.S. Balaban Correction for geometric distortion in echo planar images from B0 field variations Magn. Reson. Med., 34 (1) (1995), pp. 65-73, doi:10.1002/mrm.1910340111.
[Hutton2002](1, 2) Hutton et al., Image Distortion Correction in fMRI: A Quantitative Evaluation, NeuroImage 16(1):217-240, 2002. doi:10.1006/nimg.2001.1054.
[Huntenburg2014]Huntenburg, J. M. (2014) Evaluating Nonlinear Coregistration of BOLD EPI and T1w Images. Berlin: Master Thesis, Freie Universität. PDF.
[Treiber2016]Treiber, J. M. et al. (2016) Characterization and Correction of Geometric Distortions in 814 Diffusion Weighted Images, PLoS ONE 11(3): e0152472. doi:10.1371/journal.pone.0152472.
[Wang2017]Wang S, et al. (2017) Evaluation of Field Map and Nonlinear Registration Methods for Correction of Susceptibility Artifacts in Diffusion MRI. Front. Neuroinform. 11:17. doi:10.3389/fninf.2017.00017.