SDCFlows defines a clear :abbr:`API (application programming interface)` that divides the problem of susceptibility distortions (SD) into two stages:
- Estimation: the MRI acquisitions in the protocol for :abbr:`SD (susceptibility distortions)` are discovered and preprocessed to estimate a map of B0 non-uniformity in Hz (\Delta B_0). The theory behind these distortions is well described in the literature ([Jezzard1995], [Hutton2002]), and further discussed below (see a summary in :numref:`fig-1`). SDCFlows builds on freely-available software to implement three major strategies for estimating \Delta B_0 (Eq. \eqref{eq:fieldmap-1}). These strategies are described below, and implemented within :py:mod:`sdcflows.workflows.fit`.
- Application: the B-Spline basis coefficients used to represent the estimated \Delta B_0 map mapped into the target :abbr:`EPI (echo-planar imaging)` scan to be corrected, and a displacement field in NIfTI format that is compatible with ANTs is interpolated from the B-Spline basis. The voxel location error along the :abbr:`PE (phase-encoding)` will be proportional to \Delta B_0 \cdot T_\text{ro}, where T_\text{ro} is the total readout time of the target :abbr:`EPI (echo-planar imaging)` (:numref:`fig-1`). The implementation of these workflows is found in the submodule :py:mod:`sdcflows.workflows.apply`.
Susceptibility distortions in a nutshell
BIDS Specification
See the section Echo-planar imaging and *B0* mapping.
The displacement suffered by every voxel along the :abbr:`PE (phase-encoding)` direction can be derived from Eq. (2) of [Hutton2002]:
\Delta_\text{PE} (i, j, k) = \gamma \cdot \Delta B_0 (i, j, k) \cdot T_\text{ro},
\label{eq:fieldmap-1}\tag{1}
where \Delta_\text{PE} (i, j, k) is the voxel-shift map (VSM) along the :abbr:`PE (phase-encoding)` direction, \gamma is the gyromagnetic ratio of the H proton in Hz/T (\gamma = 42.576 \cdot 10^6 \, \text{Hz} \cdot \text{T}^\text{-1}), \Delta B_0 (i, j, k) is the fieldmap variation in T (Tesla), and T_\text{ro} is the readout time of one slice of the :abbr:`EPI (echo-planar imaging)` dataset we want to correct for distortions.
Let V represent the fieldmap in Hz (or equivalently, voxel-shift-velocity map, as Hz are equivalent to voxels/s), with V(i,j,k) = \gamma \cdot \Delta B_0 (i, j, k), then, introducing the voxel zoom along the phase-encoding direction, s_\text{PE}, we obtain the nonzero component of the associated displacements field \Delta D_\text{PE} (i, j, k) that unwarps the target :abbr:`EPI (echo-planar imaging)` dataset:
\Delta D_\text{PE} (i, j, k) = V(i, j, k) \cdot T_\text{ro} \cdot s_\text{PE}.
\label{eq:fieldmap-2}\tag{2}
BIDS Specification
See the section Types of fieldmaps - Case 3: Direct field mapping in the BIDS specification.
Some MR schemes such as :abbr:`SEI (spiral-echo imaging)` can directly reconstruct an estimate of the fieldmap in Hz, V(i,j,k). These fieldmaps are described with more detail here.
BIDS Specification
See the section Types of fieldmaps - Case 2: Two phase maps and two magnitude images in the BIDS specification.
Some scanners produce one phasediff map, where the drift between the two echos has
already been calculated, see the section
Types of fieldmaps - Case 1: Phase-difference map and at least one magnitude image.
The fieldmap variation in T, \Delta B_0 (i, j, k), that is necessary to obtain \Delta_\text{PE} (i, j, k) in Eq. \eqref{eq:fieldmap-1} can be calculated from two subsequent :abbr:`GRE (Gradient-Recalled Echo)` echoes, via Eq. (1) of [Hutton2002]:
\Delta B_0 (i, j, k) = \frac{\Delta \Theta (i, j, k)}{2\pi \cdot \gamma \, \Delta\text{TE}},
\label{eq:fieldmap-3}\tag{3}
where \Delta \Theta (i, j, k) is the phase-difference map in radians, and \Delta\text{TE} is the elapsed time between the two GRE echoes.
For simplicity, the «voxel-shift-velocity map» V(i,j,k), which we can introduce in Eq. \eqref{eq:fieldmap-2} to directly obtain the displacements field, can be obtained as:
V(i, j, k) = \frac{\Delta \Theta (i, j, k)}{2\pi \cdot \Delta\text{TE}}.
\label{eq:fieldmap-4}\tag{4}
This calculation is further complicated by the fact that \Theta_i (and therefore, \Delta \Theta) are clipped (or wrapped) within the range [0 \dotsb 2\pi ). It is necessary to find the integer number of offsets that make a region continuously smooth with its neighbors (phase-unwrapping, [Jenkinson2003]).
:abbr:`PEPOLAR (Phase Encoding POLARity)` techniques
BIDS Specification
See the section Types of fieldmaps - Case 4: Multiple phase encoded directions ("pepolar").
Alternatively, it is possible to estimate the field by exploiting the symmetry of the
distortion when the :abbr:`PE (phase-encoding)` polarity is reversed.
SDCFlows integrates two implementations based on FSL topup [Andersson2003],
and AFNI 3dQwarp [Cox1997].
By default, FSL topup will be used.
Many studies acquired (especially with legacy MRI protocols) do not have any information to estimate susceptibility-derived distortions. In the absence of data with the specific purpose of estimating the B_0 inhomogeneity map, researchers resort to nonlinear registration to an «anatomically correct» map of the same individual (normally acquired with :abbr:`T1w (T1-weighted)`, or :abbr:`T2w (T2-weighted)` sequences). One of the most prominent proposals of this approach is found in [Studholme2000].
SDCFlows includes an (experimental) procedure, based on nonlinear image registration with ANTs' symmetric normalization (SyN) technique. This workflow takes a skull-stripped :abbr:`T1w (T1-weighted)` image and a reference :abbr:`EPI (Echo-Planar Imaging)` image, and estimates a field of nonlinear displacements that accounts for susceptibility-derived distortions. To more accurately estimate the warping on typically distorted regions, this implementation uses an average B_0 mapping described in [Treiber2016]. The implementation is a variation on those developed in [Huntenburg2014] and [Wang2017].
The process is divided in two steps.
First, the two images to be aligned (anatomical and one or more :abbr:`EPI (echo-planar imaging)` sources) are prepared for
registration, including a linear pre-alignment of both, calculation of a 3D :abbr:`EPI (echo-planar imaging)` reference map,
intensity/histogram enhancement, and deobliquing (meaning, for images where the physical
coordinates axes and the data array axes are not aligned, the physical coordinates are
rotated to align with the data array).
Such a preprocessing is implemented in
:py:func:`~sdcflows.workflows.fit.syn.init_syn_preprocessing_wf`.
Second, the outputs of the preprocessing workflow are fed into
:py:func:`~sdcflows.workflows.fit.syn.init_syn_sdc_wf`,
which executes the nonlinear, SyN registration.
To aid the Mattes mutual information cost function, the registration scheme is set up
in multi-channel mode, and laplacian-filtered derivatives of both anatomical and :abbr:`EPI (echo-planar imaging)`
reference are introduced as a second registration channel.
The optimization gradients of the registration process are weighted, so that deformations
effectively possible only along the :abbr:`PE (phase-encoding)` axis.
Given that ANTs' registration framework performs on physical coordinates, it is necessary
that input images are not oblique.
The anatomical image is used as fixed image, and therefore, the registration process
estimates the transformation function from unwarped (anatomically correct) coordinates
to distorted coordinates.
If fed into antsApplyTransforms, the resulting transform will effectively unwarp a distorted
:abbr:`EPI (echo-planar imaging)` given as input into its unwarped mapping.
The estimated transform is then converted into a B_0 fieldmap in Hz, which can be
stored within the derivatives folder.
!DANGER!
This procedure is experimental, and the outcomes should be scrutinized one-by-one and used with caution. Feedback will be enthusiastically received.
There exist some alternative options to estimate the fieldmap, such as mapping the point-spread-function [Zaitsev2004], or by means of nonlinear registration of brain surfaces onto the distorted :abbr:`EPI (echo-planar imaging)` data [Esteban2016].
The workflows provided by :py:mod:`sdcflows.workflows.fit` make use of several utilities. The cornerstone of these tools is the fieldmap representation with B-Splines (:py:mod:`sdcflows.interfaces.bspline`). B-Splines are well-suited to plausibly smooth the fieldmap and provide a compact representation of the field with fewer parameters. This representation is also more accurate in the case the images that were used for estimation are not aligned with the target images to be corrected because the fieldmap is not directly interpolated in the projection, but rather, the position of the coefficients in space is updated and then the fieldmap reconstructed.
:py:mod:`sdcflows.workflows.apply` contains workflows to project fieldmaps represented by B-Spline basis into the space of the target :abbr:`EPI (echo-planar imaging)` data.
To ease the implementation of higher-level pipelines integrating :abbr:`SDC (susceptibility distortion correction)`, SDCFlows provides :py:func:`sdcflows.utils.wrangler.find_estimators`.
BIDS Specification
See the section Expressing the MR protocol intent for fieldmaps - Using *B0FieldIdentifier* metadata.
If one or more B0FieldIdentifier(s) are set within the input metadata (following BIDS' specifications),
then corresponding estimators will be built from the available input data.
BIDS Specification
See the section Expressing the MR protocol intent for fieldmaps - Using *IntendedFor* metadata.
In the case no B0FieldIdentifier(s) are defined, then SDCFlows will try to identify as many fieldmap
estimators as possible within the dataset following a set of heuristics.
Then, those estimators may be linked to target :abbr:`EPI (echo-planar imaging)` data by interpreting the
IntendedFor field if available.
If none of the two previous options yielded any workable estimation strategy, and provided that
the argument fmapless is set to True, then :py:func:`~sdcflows.utils.wrangler.find_estimators`
will attempt to find :abbr:`BOLD (blood-oxygen level-dependent)` or :abbr:`DWI (diffusion-weighted imaging)`
instances within single sessions that are consistent in :abbr:`PE (phase-encoding)` direction and
total readout time, assuming they have been acquired with the same shimming settings.
If one or more anatomical images are found, and if the search for candidate :abbr:`BOLD (blood-oxygen level-dependent)` or :abbr:`DWI (diffusion-weighted imaging)` data yields results, then corresponding fieldmap-less estimators are set up.
It is possible to force the fieldmap-less estimation by passing force_fmapless=True to the
:py:func:`~sdcflows.utils.wrangler.find_estimators` utility.
| [Jezzard1995] | Jezzard, P. & Balaban, R. S. (1995) Correction for geometric distortion in echo planar images from B0 field variations. Magn. Reson. Med. 34:65–73. doi:10.1002/mrm.1910340111. |
| [Hutton2002] | (1, 2, 3) Hutton et al., (2002) Image Distortion Correction in fMRI: A Quantitative Evaluation, NeuroImage 16(1):217-240. doi:10.1006/nimg.2001.1054. |
| [Jenkinson2003] | Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping algorithm. MRM 49(1):193-197. doi:10.1002/mrm.10354. |
| [Andersson2003] | Andersson, J. (2003) How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage 20:870–888. doi:10.1016/s1053-8119(03)00336-7. |
| [Cox1997] | Cox, R. (1997) Software tools for analysis and visualization of fMRI data. NMR Biomed. 10:171–178, doi:10.1002/(sici)1099-1492(199706/08)10:4/5%3C171::aid-nbm453%3E3.0.co;2-l. |
| [Studholme2000] | Studholme et al. (2000) Accurate alignment of functional EPI data to anatomical MRI using a physics-based distortion model, IEEE Trans Med Imag 19(11):1115-1127, 2000, doi: 10.1109/42.896788. |
| [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. |
| [Huntenburg2014] | Huntenburg, J. M. (2014) Evaluating Nonlinear Coregistration of BOLD EPI and T1w Images, Berlin: Master Thesis, Freie Universität. |
| [Zaitsev2004] | Zaitsev, M. (2004) Point spread function mapping with parallel imaging techniques and high acceleration factors: Fast, robust, and flexible method for echo-planar imaging distortion correction, MRM 52(5):1156-1166. doi:10.1002/mrm.20261. |
| [Esteban2016] | Esteban, O. (2016) Surface-driven registration method for the structure-informed segmentation of diffusion MR images. NeuroImage 139:450-461. doi:10.1016/j.neuroimage.2016.05.011. |