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}}$